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ABSTRACT 


Three methods for solving the dynamical equations of flutter 
are presented. The application of these procedures to systems 
described by more than three degrees of freedom is discussed. 
It may be pointed out in this connection that the number of 
degrees of freedom is not always a true indication of the complex- 
ity of a problem. Problems are more legitimately compared 
when the dynamical equations are reduced to a standard ‘‘canoni- 
cal’’ form. 

The essential steps in the methods presented are illustrated 
with detailed numerical examples. Moreover, a knowledge of 
the mathematics required in the theory is not presumed. All 
such material is developed in a preliminary section of the paper. 

Convenient formulas for solving the general cubic and quartic 
equations are given. The formulas have been found to be well 
suited to machine computation. A concise tabular process for 
solving a system of linear equations is explained with a numerical 
example. 

Practical suggestions, based on experience, are offered relative 
to the organization of flutter calculations. The results of an ex- 
tended time study are presented which compare the methods on 
the basis of computing man-hours required. This comparison 
is rounded out with a discussion of the relative merits of the 
methods with respect to other considerations that are not in- 
cluded in the time study. On the basis of these results defi- 
nite recommendations concerning the use of the methgds are 
made. 


INTRODUCTION 


Oz OF THE MOST TROUBLESOME PROBLEMS confront- 
ing the flutter analyst is that of expeditiously 
solving the dynamical equations once they have been 
obtained. This problem is daily becoming more im- 
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portant because of the increased use of analyses involv- 
ing four, five, six, and even more degrees of freedom. 
With these thoughts in mind, an investigation of the 
various analytical methods that seemed applicable to 
the flutter problem was undertaken. The most per- 
tinent results of that study are presented in this 
paper. 

It is evident that this problem could also be studied 
from other points of view. Thus the possibility of 
constructing an equivalent electrical circuit could be 
investigated. The use of “punched cards’’ would also 
be worthy of consideration. However, a preliminary 
survey of the analytical methods of solution seemed to 
be a logical prerequisite to any other type of research. 


MATHEMATICAL PRELIMINARIES 


It has been thought desirable to give a preliminary 
survey of the mathematical conventions, definitions, 
and theorems to be used throughout the paper. The 
material to be presented covers the essentials of deter- 
minants and linear equations. 

The basic notation of tensor analysis is used in this 
paper. This notation makes the use of matrices super- 
fluous. Nevertheless, it has been thought convenient 
to indicate the correlation with the usual matrix ter- 
minology from time to time. However, all remarks 
pertaining to matrices may be disregarded without im- 
pairing the continuity of the development. 


Summation Convention 


Whenever an index occurs twice in a single term, it 
shall be understood-that the term represents the sum 
of the individual terms formed by replacing the repeated 
index with the positive integers 1 to n. 


Diya; = Dua, + Did, +...+ Dindn (A) 
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It is to be noted that always represents a definite 
positive integer when used in this connection; it is 
never to be used as an index calling for summation. 

It follows as a direct consequence of the summation 
convention that 


= Dum =.. = Duta (1) 


Eqs. (1) state that the repeated index used is imma- 
terial, provided, of course, that it does not conflict with 
the use of the same index in a different sense. For this 
reason, j, k, a as used in Eqs. (1) are termed “dummy” 
indexes. 

Consider now the equality 


Diya; = = 23 (2) 


It is to be noted that this equality represents a set of n 
equations; the ith equation of the set may be written 
as 


+ + + Dindn = + +... 
+ Zbindn (3) 


It should be observed that the same ‘“‘free’”’ index 7 oc- 
curs in both members of this equation. This illustrates 
the general rule that in any equation all terms should 
agree as to free indexes. This important check should 
always be made when manipulating equations. 

A second convention exists in the statement that all 
free indexes shall be understood to represent the posi- 
tive integers from 1 to n. 

In certain expressions it is convenient to use a re- 
peated index more than twice. The summation con- 
vention will be extended to cover such cases; thus it 
is to be understood that 


On the other hand, it is sometimes convenient to waive 
the summation convention with respect to an index. 
An explicit statement will be made in all such cases; 
for example, the following equations will occur: 


D = (5) 


not summed on a. It may be observed that there is 
actually no ambiguity in this example since both 7 and 
a are free indexes in the left member of the equation 
and must therefore be regarded as free indexes in the 
right member. 

The function 6,,, known as ‘“Kronecker’s delta,’’ will 
be used extensively; it is defined by the equations 


6, = 0, fort (6) 
6, = 1, fort = 7. (7) 


In matrix terminology 5,; would be the identity matrix 
I. 

Let Eqs. (2) be reconsidered in the light of Eqs. (6) 
and (7); carrying out the summation on j in the right- 
hand member of Eqs. (2) yields the result 


Digay = Za; (8) 


Linear Equations and Determinants 


“The general system of m linear equations in the n 
unknowns x, may be written as 


= hy (9) 


If all the h; are zero, the equations are called “‘homo- 
geneous”; otherwise the equations are “‘nonhomogene- 
ous.” The typical set of equations for = 3 is 


+ Zire + gists = 
+ + = he (10) 


+ + = hs 


Definition 1 

A determinant is a square array of numbers to which 
a single number, called the value of the determinant, is 
assigned according to a rule to be stated subsequently. 
The individual numbers of the array are the elements, 
the horizontal sets of elements are the rows, and the 
vertical sets of elements are the columns of the deter- 
minant. The number of rows (or columns) is the order 
of the determinant. The diagonal set of elements ex- 
tending from the upper left-hand corner to the lower 
right-hand corner is the main diagonal of the determi- 
nant. The typical element of the determinant is written 
as 21, 7 referring to the row and j referring to the column 


containing the element. The determinant itself, as_ 


well as its numerical value, is designated by the symbol 


The coefficients of Eqs. (10) form a determinant: 


gi gis 
£21 S22 
Su 


(B) 


The elements (gu, Z12, g:3) form its first row; note that 
this row may be denoted by gi;. Similarly, the second 
column is gz, etc. : 

The numerical value of a determinant depends upon 
certain related quantities GY, which are now defined. 


Definition 2 


The cofactor G” of the element g,; is the determinant 
formed from |g;,| by crossing out the ith row and jth 
column of |g,,| and giving the resulting determinant 
of order (n — 1) the sign (—1)**?, 

Considering determinant (B) it is found that 


= (-1)!# 823 

(C) 
G8 = (—1)2+3/8 

£31 £32 


The numerical value of a determinant may be found 
by applying the following theorem. 
Theorem 1 
leul = GY gy 
(summed on either 7 or j7, but not both). 


it 


it 


METHODS FOR SOLVING THE FLUTTER EQUATIONS 5 


Stated verbally, the value of a determinant may be 2 db). 
found by each element of a row (or column) values of determinants of the type 
by its cofactor and then summing the resulting products. given indicates that 

Of course, since each cofactor is itself a determinant 
of order (n — 1), the same theorem applies to it. Thus 
the evaluation of less finally rests on the numerical 


a b 


= od be 


Let determinant (B) be expanded by Theorem 1; thus, using the first column of |g,,|, it follows that 


ges] = gu(geages — — — + — Lissa) 
Consider now a determinant of the type 


Qu + by 
A= + bay Qe 
+ bs 


Expanding in terms of the elements of the first column it is evident that 


bu ay ay 
A= Ax + 
bs 


Applying this result repeatedly, it is possible to express an mth order determinant, each of whose elements is the 
sum of two quantities, as the sum of 2" determinants of the mth order in each of which the quantities occur 
singly. 

It may be verified by direct expansion that 


u + du nu dy u dy bu bu 


The result for the general case may be put into a neat and easily remembered form. As an example, consider 
the case m = 4 for which the determinant is 


dy be ais t+ dis au + 

+ da + dae + + Dag 

Let two auxiliary determinants be defined: 

Qu Gy ay bu bis by 

A. = Gog and Ay = dae bas da 


Moreover let A,(ij ... k) be the determinant formed from A, by replacing the columns i, j, ..., k of A, by the 
corresponding columns 7, j7, ..., R of A,; an analogous definition holds for A,(ij ... k). Using these definitions 


‘it may be shown that 


A= Ad(1) + Ad(2) + + Ad(4) + Ac(12) + A,(18) + A,(14) + AQ(23) + AQ(24) + A.(34) + 
A4(123) + Ag(124) + A,(134) + A,(234) + A,(1,234) (12) 


It will be noted that A,(1,234) = A,; similarly it is seen that A,(123) = ,(4); etc., so that the expression 


. given for A may be modified if desired. The extension of the expression given for A to the general case is straight- 


forward but will not be given. 
The following theorem states a fundamental property of determinants. 


5 
) 
or 
h 
s 
e 
¢ 
r 
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Theorem 2 


GY gx; = 0, for ixk 
G%g, = 0, forj #k 


Thus the sum of the products of the elements of a 
row (or column) by the cofactors of the corresponding 
elements of another row (or column) is zero. 

With respect to determinant (B) Theorem 2 states 
that 


+ guG*? + = 0 (E) 
This fact can be verified by direct substitution; thus 
GY = (—1)1+2/8% 
833 
G2 = (—1)2+2/82 
£31 £33 
G? = 
£2 823 


so that the foregoing sum becomes a 
—gu(g8ss — ages) + gei(gugss — — 


Theorems I and 2 can be conveniently summarized 
with the aid of the following definition. 


Definition 3 
g” = cofactor of in the determinant|g,,| 
lg ul 
The quantities g“ are called the ‘‘normalized’”’ cofactors 
of the elements g,;;. The g are defined only when 
\gi| * 0. Theorems 1 and 2 are now expressed by the 
equations 


= On (13) 
= (14) 


These results may be used to obtain the solution of a 
system of linear nonhomogeneous equations; thus let 
Eqs. (9) be rewritten as 


= hy (15) 
It is now assumed that 
(16) 


Multiplying Eqs. (15) by g* and using Eqs. (14), it is 
found that 


= (17) 
Summing the left members of Eqs. (17), it follows that 
= (18) 


_ The solution of Eqs. (15) given by Eqs. (18) is some- 
times known as Cramer’s rule. 
Consider now a set of homogeneous equations: 


= 0 (19) 
Egs. (19) clearly have the “‘trivial” solution x; = 0 for 


j = 1, 2, ..., m regardless of the nature of the deter- 
minant gus. This is the solution given by Eqs. (18) 
on the assumption that lgus| ~ 0. If, however, 


=0 (20) 


it will be shown that Eqs. (19) possess nontrivial solu- 

tions in addition to the trivial solution mentioned. 

This result is of fundamental importance for all vibra- 

tion work. The proof that follows is rather lengthy 

but has been included because the necessary theory 

has already been developed in the preceding pages. 
Let Eqs. (19) be partially expanded as follows: 


+ Lintn = =1,2,...,.8—1 (21) 
+ =O) =1,2,...,m —1 (22) 

Eqs. (21) may now be rewritten as 
= = 1,2,...,m —1 (23) 


It should be noted that the determinant formed by the 
coefficients of the x,;/x, in Eqs. (23) is the cofactor G" 
of the element ga» in the determinant |g,,|._ The ratios 
x,/X, may now be determined from Eqs. (23) by Cra- 
mer’s rule on the assumption that 


G™ #0 (24) 

Let y* represent the normalized cofactors of the deter- 

minant G™; applying now the solution given by Eqs. 
(18) to Eqs. (23) it is found that 

= 4,7 = 1,2,...,m —1 (25) 


The quantities x,/x, are, of course, solutions of Eqs. 
(21); in order that they be solutions of Eqs. (19), it 
is necessary that they satisfy Eqs. (22) as well. Sub- 
stituting from Eqs. (25) into Eq. (22), it follows that 


— + Sante, = 0; 4 = 1,2, ...,8 — 1 

(26) 
Since it has already been assumed that x, ¥ 0, Eq. 
(26) can be true if and only if 


Zin + = 0; 1, j =1, 2, 1 (27) 


It will now be shown that Eq. (27) is a logical conse- 
quence of Eq. (20) together with Theorems 1 and 2. 
Thus from Theorem | it may be stated that 


les] = G%e-s (not summed on j) (28) 
By virtue of Eq. (20) this becomes 
G"g,; = 0 (not summed on /) (29) 
By Theorem 2 it is known that 
= 0, forj k (30) 


Combining Eqs. (29). and (30), therefore, it may be 


stated, without any restriction on the free indexes, that 


= 0 (31) 
In particular it may be stated that 


(31) 
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gn = (32) 


In Eqs. (28) through (32) it is understood that all 
indexes range from 1 to m; it is now convenient to re- 
write Eqs. (32) as 


Gg + = 0 (33) 


In Eqs. (33) the summation index s ranges from 1 to 
n — 1 while the free index k goes from 1 tom. From 
Egs. (33) it follows that : 


gu = $= 1,2, (34) 


Substituting for —g,, into Eq. (27) from Egs. (34) it is 
seen that 


(G'"g55/G"") "gin + = 0 (35) 


It is to be noted that each of the indexes 7, j, s ranges 
from 1 to m — 1 in Eq. (35). It is now recalled that 


= by; (36) 
Hence Eq. (35) becomes 


gin + = 0 (37) 
Summing on s, it is now seen that 
(G"/G"")gin + = 0 (38) 
Hence 
G"gin + SnnG™" = 0 (39) 


However, this last statement is correct since the left 
numberof the equation is | Zul, which is zero. 

It therefore follows that Eq. (27) is valid; thus 
Eqs. (25) constitute the solution of Eqs. (19). 


GENERAL THEORY 


The basic dynamical equations of flutter have been 
derived by several investigators. As already stated, 
this paper is concerned primarily with the problem of 
solving the equations expeditiously. It may be stated, 
however, that the general approach to the flutter prob- 
lem itself is that of Bleakney. In the notation of 
Bleakney, the dynamical equations to be considered 
are: 


[(m Aw) + ies) 0 (1) 


(not summed on 7). 
The indexes i and j each range from 1 tom. The j 


preceding the g,,; represents the V—-1. For conven- 
ience the following substitutions are made: 


— Ay=Cy; miyQy = ky (2) 


Using the definitions given by Eqs. (2) it is seen that 
Eqs. (1) become 


{Cy [(wo?/ w?)(1 + 5813) la, =0 (3) 


In theory the damping parameters g;; are assumed to 
be known numbers. Actually, they are, in general, 
unknown and must be experimentally determined. An 
alternative approach presents the flutter stability in 
the form of stability curves exhibiting the “equivalent 
damping” factor g required for critical stability of the 
system as a function of air speed. This is mathe- 
matically equivalent to putting 


gy =g for all i and j (4) 
and then introducing the complex variable 


Z = (wo?/w*)(1 + jg) (5) 


into Eqs. (3). Making the substitutions just given, it 
is found that Eqs. (3) become 


(Ciy — = 0 (6) 
Assuming now that 
(7) 


it is possible to multiply Eqs. (6) by the normalized 
cofactors thus 


(k™C,, — Zk*k;,;)a; = 0 (8) 
Performing the summation on 7 and defining 
Day = (9) 
yields the result 
(Das — Zbqs)a3 = 0; a,f7 = 1,2,...,m (10) 


Eqs. (10) will be referred to as the “canonical” form 
of the flutter equations. 

Before proceeding with the solution of Eqs. (10) it is 
desirable to consider the reduction of other systems of 
equations which often arise to the canonical form just 
given. For example, when rigid body motions of the 
airplane are included in a flutter analysis as separate 
degrees of freedom, it becomes necessary to consider a 
system of equations 


(Cy Zkij)a; = 0 (11) 
where now 
= 0 (12) 


due to the fact that m equations, say, of the set are of 
the form 


=0 (13) 


It will be noted that, since |k,,| = 0, the reduction to 
canonical form previously given is not valid. How- 
ever for this case (n — m) of the a, may be algebraically 
eliminated to bring the equations into a form for which 
the previous reduction is valid. This preliminary elimi- 
nation of variables will be explicitly carried through 
for the cases (n — m) = 1 and (n — m) = 2. For the 
case (n — m) = 1 it is convenient to rewrite Eqs. (11) 
and (13) in the forms 


7 
ter- 
(18) 
solu- 
ned. 
ibra- 
gthy 
eory 
cox) 
(22) 
| 
y the 
4 G™, 
Cra- 
(24) 
eter- 
Eqs. 
(25) 
9), it 
Sub- 
that 

-1 
(26) 
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(30) 
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(Cy = = = (14) 
= 0 j=1,2,....%*% +1 (15) 
From Eq. (15) it follows that 
nti 


Thus writing Eqs. (14) in the form 
(Cia — Zhig)a + (Cint1 — = 0 (17) 


and noting that kin+1 = 0, since ky+1, = 0, it follows 
by direct substitution that 


nt+1 
1,2, ...,” (18) 
Defining 


Big = — (Contr Cotra/ Cots (19) 
it is clear that 
(By — = 0; =1,2,...,n (20) 


Evidently the reduction to canonical form is now 
identical to that used with Eqs. (6). 14 will be noted 
that, once the a, of Eqs. (20) have been determined, 
@,+1 can be found from Eq. (16). 

For the case (n — m) = 2 the equations are 


(Cry — = (21) 

(22) 

From Eqs. (22) it is found that 

and 


Substitution of a,+2 from Eq. (24) into Eq. (23) yields 
the result 
— Cat2 nt2Catta (25) 
Defining the quantity in parentheses to be tia, Eq. (25) 
becomes 


Onti = =1,2,...,” (26) 
‘Substituting for a,+1 from Eq. (26), Eq. (24) yields 
the result } 
( Cate ntiCatia — Cat ntiCat20 (27) 
nti Cate nt2 — Cati ati 


Thus 


Ont2 = =1,2,...,0 (28) 


where f., represents the quantity in parentheses. 


Egs. (21) may be expanded to read - 


(Cia — + (Cine: — + 
(Cint2 — ZRint2)@xr2 = 0 (29) 


Since = 0 and = 0, it foHlows, from: the 
fact that ky = ky, that 


Rint: =0; = 0 (30) 


Substituting from Eqs. (26), (28), amd (30) imto 
Eqs. (29), it is found that 


{(Cia + traCints + taCint2) — = 0 (31) 
Making the definitions 
By = Cia + + 4, a = 1,2, ..., 
(32) 
it ts evident that Eqs. (31) become 
(Bi, — = 0; 4,2 =1,2,...," (33) 


The reduction to the canonical form of Eqs. (10) is 
now evident. Once the a, of Eqs. (33) have been de- 
termined, and @,+2 are found from Eqs. (26) and 


(28), respectively. 
It should be noted that since ky, = ky, it follows that 
¢ =k (34) 
It may also be recalled that 


where the K® are the cofactors of the k,,. 

Before proceeding with the general theory it seems 
desirable to give a numerical illustration of the reduc- 
tion of Eqs. (6), say, to canonical form. 


Illustration of Reduction of Equations to Canonical Form 
Consider the equations 


(4 — 2Z)a, + (5 — + = 
(1 — + (3 + + = 0> (1) 
5a, + 2a2 + (7 — 3Z)as = 0 
Evidently the k,, are 
2 1 0 
ky => 1 —1 0 (2) 
.4 8 
Recalling the definition of cofactor it will be seen that 
K¥%=-3 K*®=-3 K*®= 0 3 
8) 
Also it is clear that , 
= = —9 (4) 


Hence, using Eqs. (34) and (35) yields the result 


1/, —2/, = (5) 
M= 0 0 
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Substituting into Eqs. (9) gives the values 


Du = 1/3(4) + + 0(5) = 
Dy = 1/3(5) + 1/2(3) + 0(2)= 
Dig = 1/3(6) + 1/3(4) + 0(7) = 
Du = 1/3(4) — + = 
De = 1/3(5) — */x(3) + 0(2) = (6) 
Dy = 1/3(6) — + 0(7) = —*/s 
Dun = 0(4)+ O(1) +1/:(5) = 
Dx = 0(5) + 0(3) +1/3(2) = */s 
Ds = 0(6) + 0(4) +1/:(7) = "/sj 


Thus Eqs. (10) become 


— Z)ai + 8/sa2 + = 0 
+ (—1/3 — Z)d2 — = 0 (7) 


The solution of Eqs. (1) has thus been reduced 
to the solution of the canonical equations given 
above. 

It has been shown that in all cases the dynamical 
equations can be reduced to the canonical form typified 
by Eqs. (10). In all that follows, it is understood 
that the equations are in canonical form. For con- 
veniertce of reference, Eqs. (10) are rewritten here 
as 


(Dy = 0; 4,j =], (36) 


These equations are linear and homogeneous in the 
unknowns a;; hence solutions other than the trivial 
solution a, = 0 will exist if and only if the determinant 
of the coefficients is zero—i.e., 


— = 0 (37) 


Eq. (37) is commonly known as the ‘frequency equa- 
tion.” It is also called the ‘‘characteristic equation”’ 
of the matrix D,,._ For flutter work the name ‘“‘stabil- 
ity equation” is also used. The determinant in Eq. 
(37) is known as the “‘stability determinant’; in the 
classical literature it is the ‘‘characteristic determinant” 
of the dynamical matrix D,;. 
gree m in Z obtained by expanding the stability deter- 
minant will be called the “stability polynomial.” ‘The 
roots Z, of the stability equation have various designa- 
tions: ‘‘latent roots” of the matrix D,,; ‘‘eigenvalues’’ 
or ‘‘characteristic values’’ of the system. The set of 
a, obtained from Eqs. (36) when Z, is substituted for Z 
will be denoted by a,,; the quantities a,, constitute 
what is called the “eigenfunction” or ‘characteristic 
function’”’ belonging to the eigenvalue Z,. 

In matrix language, the a, constitute the “modal 
column” corresponding to the “latent root” Z,. For 
flutter work it is convenient to label the a,, the ‘“‘ampli- 
tude ratios’”’ associated with Z,. This emphasizes the 
fact that only the ratios of the unknowns ay, are de- 
terminable from Eqs. (36), since the equations are 
homogeneous. 


The polynomial of de- . 


So much now for notation and terminology; the 
method of iteration will now be developed. The first 
step will be to establish certain relationships among 
the a,;,. For this purpose consider the associated 
equations 


— 54Z) = 0 (38) 
[Dy — = 0 (39) 


Remembering that a determinant is unchanged in 
value by the interchange of corresponding rows and 
columns, it becomes clear that Eq. (39) is identical to 
Eq. (37). However, it is to be noted that since, in gen- 
eral 


Diy Dy (40) 


Eqs. (36) and (38) are not identical sets of equations. 
The positions of the a; and A, in Eqs. (36) and 
(38) correspond to the observation that the a, con- 
stitute a column matrix, while the A; form a row ma- 
trix. 

Let Ag, denote the solution of Eqs. (38) correspond- 
ing to the eigenvalue Z,;—i.e., 


— = 0 (41) 
Referring now to Eqs. (36) it is known that 
(Diy — = 0 (42) 


From Eqs. (41) and (42) it follows that 


AgyDy, = Z,Ag; (not summed on 
= Z Aig (not summed on a) 


Multiply the first of these equations by a,, and the 
second by Ag;; thus 


ji, = (not summed on 8) (43) 
= (not summed on a) (44) 


Observe that since both 7 and j are dummy indexes, the 
left members of Eqs. (43) and (44) are identical; there- 
fore 


(Za — Zs)Agdig = 0 (45) 

Assuming now that 
Z, # ita B (46) 
it follows that 
Agidig = 0, for a ¥ B (47) 


This is the so-called orthogonality property of the eigen- 
functions of Eqs. (36) and (38). Using Eq. (47) it 
follows from either of Eqs. (43) or (44) that 


AgD ij. = 0, for a B (48) 


Let it now be assumed, that the Ay, say, have been 
determined in some way; using Eq. (47), the number of 


the 
(30) 
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(33) 
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equations to be solved can now be reduced from to 
nm — 1. The method of reduction has already been il- 
lustrated with Eqs. (14) and (15). The process will, 
however, be repeated at this point. 

For 8 = 1, Eq. (47) may be written as 


Aja; = 0 (49) 


with the understanding that the a; corréspond to 
eigenvalues other than Z;. From Eq. (49) it follows 
that 


a, = = 1,2,...,m —1 (50) 

Eqs. (36) may be partially expanded to read 
(Dis — + (Din — = 0 (51) 
— + (Dan — = 0 (52) 


where i,j = 1,2, ...,m — 1. Substituting into Eqs. 


(51) from Eq. (50) and noting that 6,, = 0, it is found 


that 


[ (2. fa, = 0; 4,j =1, 2, — 1 
In 
(53) 


From the manner in which Eqs: (53) follow from Eqs. 
(51) it is clear that the roots of 


|[Dis — (A1j/Ain)Dinl — = 0 (54) 


are the roots of Eq. (37) with the exception of Z;. The 
dj, corresponding to the roots Z, of Eq. (54) may be ob- 
tained from Eqs. (53). These values, together with a, 
determined from Eq. (50), constitute the solutions of 
Eqs. (36) corresponding to all Z, except Z; again. It 
may also be noted that Eq. (52) is redundant for the 
purpose at hand. Defining 


Diy) = Dy — (55) 
it is evident that Eqs. (53) become 
(Dij’ — 6yZ)a; = 0; 4,7 =1,2,...,m—1 (56) 


The equations to be associated with Eqs. (56) are evi- 
dently 
Aj (Dy' — = 0 (57) 


Note, however, that the A,’ of Eqs. (57) are not 
identical to the A, of Eqs. (38). ‘ 

Clearly, all that has been demonstrated with respect 
to Eqs. (36) and (38) could now be repeated using Eqs. 
(56) and (57). Thus it appears that the complete de- 
termination of the Z, and a,, corresponding to Eqs. 
(36) rests essentially on the possibility of developing a 
method for obtaining Z;, and Ai,;. These quantities 
can be found by an iterative process to be described at 


this point. The validity of the procedure will be 
demonstrated subsequently. 
Iterative Solution of the Dynamical Equations 
Consider Eqs. (36) written in the form 

= Za; (58) 
The multiplication-summation process given by the 
left member of this equation is used to construct a 
sequence of “iterates” Thus, beginning with an 


arbitrarily chosen set of m numbers a,, the sequence 
is generated as follows: 


Da; 

Da ™ = 
i (59) 

= a, 


In matrix terms this process is described as repeated 
postmultiplication of the matrix D whose elements are 
the D,, by an arbitrary column matrix g whose ele- 
ments are the a;(°). 

If the procedure defined by Eqs. (59) were to be car- 
ried out for a “‘suitable’’ numerical example, it would 
be noticed that the ratios of the elements of the rth 
iterate to the corresponding elements of the (r — 1)th 
iterate tend to a common limiting value as r increases. 
It will be clear upon consideration of Eqs. (58) that this 
common ratio, together with the last iterate used in 
obtaining it, constitutes a solution of those equations. 
It will be proved, in fact, that the common ratio may be 
taken for Z,, so that the last iterate is necessarily ay. 
Evidently the calculations just described could also be 
carried out with reference to Eqs. (38). Thus writing 


Eqs. (38) as 
ZA, (60) 


a sequence of iterates A, could be constructed as 
before, beginning with an arbitrary set of m numbers 
A,,; viz. 


= Ai”? (61) 


This corresponds to the repeated premultiplication of 
the matrix D, whose elements are the D,,, by an arbi- 
trary row matrix Vo with elements A,. The ratios 
of corresponding elements of A; and A,‘-” will 
again tend to a common value as ¢ increases provided 
certain conditions are met. As before it may be shown 
that this ratio is actually Z,, so that from Eqs. (60) 
it follows that the final iterate may be taken as 


* Ay. 


Before demonstrating the convergence of the itera- 
tion process, it has been considered desirable to illus- 
trate the procedure with a detailed numerical ex- 
ample. 
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Example of the Iterative Method 
Let the solutions of the following equations be required: 


(3 — Z)a — a: + a3 — a=0 
2a, + (8 — Z)a, — a3 — 2a, = 0 (1) 
—a + 2a, + (—1 — Z)as — a=0 
a — a: + a3 + (1 — Z)a, = 0 


The equations given are in the form of Eqs. (36); the iterative determination of Z; will be carried out along the 
lines of Eqs. (61). It is not, however, necessary to explicitly write out Eqs. (38) for this purpose. This state-’ 
ment will be clear from the work to follow; let A, be taken as the set of numbers (0, 0, 0, 1); 


A,® = 0; A, = 0; A; = 0; A =] (2) 
The substitution of A, into the first of Eqs. (61) can be conveniently tabulated as follows: 


1 =1 0 0 
0} 2 -i -2 0 0 
2-1 0 0 
1 1 =I 


oO 


1 -1 1 1 ~ AM 


It will be noted that the set of D,, is copied from the given equations; the A, are arranged vertically after 
which each of the A; is multiplied into the horizontal set of numbers to its right. The resulting products are 
then added vertically. The process is now repeated; thus 


1/ 3 1 -1 3 1 
2 8 -1 || -2 -8 1 2 
2-1 |} -1 2-1 —-1 
1} 1 —1 1 1} 1 1 1 


1 —8 2 


3 -1 1 -1 3 -1 1 -1 
-s| 2 8 -1 -64 8 16 
2 -1 -2 4 -2 -2 
i} 1-1 1 1 


-@ 8 14~ A,® 


—14/ 3 -1 1 -1 —42 14 -14 14 
2 8 -1  ||-124 -—496 62 124 
si-1° 2 || -8 14 -8 -8 
1-1 #1 1 14-14 1 14 


—160 —480 54 144 ~ A," 


3 —1 — 480 160 —160 160 
—480)| 2 8 —1 ||—-960 —3,840 480 960 
2-1 || —54 108 —54 
144 J] —-1 1 1 144 —144 144 144 


—1,350 -3,716 410 1,210 ~ 


As the process continues, attention should be given to the successive ratios A; /A;"-"; using A; and A,‘ 
it is found that 


A A 


A,® 


7.74; = 7.58; = 8.40 


It appears, therefore, that Z, has a value in the neighborhood of 8; continuing the process, A;® is found: 
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—1,350) 3 1 —1 — 4,050 1,350 —1,350 1,3 
—3,716)| 2 8 -1 [7,432 .-—20,728 3,716 7,432 
410)|—1 2-1 —410 820 —410 —410 
1,21 1 1 1 1,210 -—1,210 1,210 1,210 
—10,682 —28,768 3,166 9,582 ~ A,® 
The ratios A,/A,® are now (7.92, 7.76, 7.73, 7.98). 

Continuing once more it is found that A,” is (—83,166) 
(—83,166, —222,712, 24,502, 74,634); moreover the 
ratios are fairly constant; thus fue 

A Ag? Ay”? In the same manner it will be found that 
Dy! = —3.99 Ds! = 1.328 
Therefore Z, will be taken to be 7.77 and Ay, will be Dal = 222 Day’ = 
Dz’ = —0.99 = —0.672 


taken to be A;. Using A1,, a, will be eliminated from 
the given equations as described previously; compare 
Eqs. (50), (55), and (56). 
From Eggs. (55) it is known that 
Diy! = Dy — (A1j/Au)Du 


Substituting numerical values into these equations it 
follows that 


Dy’ = (— 83,166) (-1) = 38 — 1.115 = 
— 1.885 

Da’ = 2 — 2.230 = 
74,634 —0.230 


Hence the reduced equations corresponding to Eqs. 
(56) are 


(1.885 — — 3.99a2 + 1.3284; = 0 
—0.230a; + (2.02 — — 0.34403; = 0 (3) 
—2.115a; — 0.994, + (—0.672 — Z)as = 0 


Clearly, the process already illustrated can be applied 
to this reduced set of equations, thereby obtaining Z 
and Ao’. 

On the other hand, either of the two methods to 
be given subsequently may be used to complete the 
solution. 


Should the amplitude ratios corresponding to Z; be desired, the iterative procedure outlined by Eqs. (59) may 
be applied. This process can also be cast into a convenient tabular form as follows: choosing a, to be (0, 1, 
0, 0), the calculation of a," according to the first of Eqs. (59) may be arranged as shown below: 


0} 3 2 1 
8 2 
1 


8 2-1 
0 
—1 8 2 -lwa™ 


It should be noted that the rows and columns of the array used in computing the A,‘ have been interchanged; 


using this new array the process is identical to that shown for the A,”. 
large numbers in the successive iterates, the process originally used will be modified; 


However, to obviate the appearance of 
thus each number of a,‘ 


will be divided by a," which is —1. Symbolically then a,;9 = —1(1, —8, —2,1); the numbers in the paren- 
thesis will now be used to continue the iteration process. The justification of this procedure is apparent when 
it is remembered that (1, —8, —2, 1) may be regarded as a new arbitrary first approximation, a;%. Hence, 
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8 2 


-2 1 


3 —!1 1 


I 
e 
g 
t 
is 
ti 
fe 
Ww 
ni 
E 
Ww 
| of 
kr 
so 
mi 
tic 
fai 
of 
be 
pre 
tio 
| 8 —64 —16.. 8 
-2 -1 1 Mt 
foll 
8 —62 8 


2.115 


nged; 
nce of 
of a, 
paren- 


METHODS FOR SOLVING THE FLUTTER EQUATIONS 13 


The set (8, —62, —16, 8) is now written as 8(1, —7.75, —2, 1); continuing as before—i.e., regarding (1, —7.75, 


—2, 1) asa new a; — 


i]s: 2-1 1 3 2 -1 1 || 
—7.75||—1 8 2 || 7.75 -62 -15.5 7.75 
1-1 —1 2 2 

7.75 —60 —15.5 7.75 
Again (7.75, —60, —15.5, 7.75) is written as 7.75 (1, —7.74, —2, 1); thus 
3 2 -1 1 3 2 1 
8 2 7.74 —61.9 -15.48 7.74 
1 —-1 1 —2 2 2 —2 
-2 -1 —2 1 
7.74 —59.9 -—15.48 7.74 


Since the results of the last two iterations are practically identical, it is found that 
Zi = 7.74 and aa = (1, —7.74, —2, 1) 


It will be noted that Z, as determined by the A,” 
was taken to be 7.77; hence the question of accuracy 
must be considered. Certainly, for engineering pur- 
poses, either 7.77 or 7.74 is acceptable. One might 
even argue that 8 is good enough. However, such ar- 
guments ignore the essentially mathematical nature of 
the question. The important point to be remembered 


is that the A,” associated with Z; are to be used in ob- — 


taining the D,,;’ which determine the reduced Eqs. (3). 
Clearly then, the A,‘ taken as A;,; must be sufficiently 
well determined so that the D;,;’ are themselves defi- 
nitely determinate. 

In the example given, the calculations leading to 
Eqs. (3) indicate that no significant change in the D,,’ 
would result from a continuation of the iterative process 
beyond A,”. 

A few words may be in order concerning the choice 
of A, and a,;. Evidently, if Z; were approximately 
known, A, and a, would be chosen so that at least 
some of the ratios A;‘?/A, and a,/a, approxi- 
mate the known value. A consideration of the equa- 
tions to be solved sometimes makes it possible to guess 
fair approximations to the Z,. Thus if the determinant 
of the coefficients of the equations has its largest num- 
bers along the main diagonal, these numbers will ob- 
viously be fair approximations to the Z,. 


Convergence of the Iterates 


The convergence of the iteration process will now be 
proved with reference to Eqs. (36). From the defini- 
tions given for Z, and a,, it is known that 


= (not summed on a) (62) 


Multiplying Eqs. (62) by a** and summing on a@ the 
following equations may be obtained: 


D ya" = a" Z Lia 


= 
Da = a? Z 


Replacing the free index k by j for convenience, it is 
clear that 


Dy = dia (63) 
Let this expression for D,; be used in Eqs. (59); thus 
af) = (a*Z,a 
so that 
a = (a (64) 


The change of dummy indexes is made to avoid any 
possible confusion in the work to follow. Manipulat- 
ing the right member of Eq. (64), it is found that 


a; = 
a; te” Z gay 
= Zraga”Z,a, 


Thus it is found that 
= aa, (65) 
In a similar fashion it follows that 


a,® = Dyas? = (Zin) 


Thus it is seen that 
By the process of mathematical induction it may be 
shown that in the general case 
a,” = r=0,1,2,... (67) 


| = 
yplied 
ng Ze 
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Using Eqs. (67) it is clear that 


a,” /a,*-» = /Z (68) 
Tn partially expanded form this equation reads 
a” Zy'aya"!a, + Za aaa, +... +Z,'a; (69) 
+ + ... + a, na*"a, 
Factoring a Z," from the numerator and a Z,"~! from the denominator, Egs. (69) become 
r r 
+ (2) a gaa, +...+ (7) a; 
(70) 
a (r—)) r—l 
i + (2) (2) Qin 
Zi Z 1 
If now it is assumed that subject to the inequality previously given; viz.: 
\Z,/Z:| <1, fora >1 (71) |Z./Z:| <1, fora >1 (78) 
it follows from Eqs. (70) that As indicated previously, the sequence A,” has the 
limiting value Ay. Knowing Eqs. (56) can be 
= 
Him constructed. A repetition of the process used in deter- 


mo 


As already stated previously, a,” may be taken as dy. 

The convergence of the procedure outlined by Egs. 
(61) can be carried out in an analogous manner. Thus, 
since by definition 


A = ZA = ZgA gi (not summed on (73) 
it follows that 
ji = Zp A*A 
ji = Z,A**A gi 
Dui Zp A**A 
Hence it has been shown that 
Dy = Z,A%As: (74) 
Substituting this result into Eqs. (61) yields the follow- 
ing equations: 


A® Z,A%A As 
Ai = ZA 
Through successive substitutions it is found that in 
the general case 


A, = 


(75) 


(76) 
As before it will be seen that 
lim (A,®/A,“-») = Z,;; 4 = 1,2, ...," 


mining Z, and A, will now yield Z; and A2,’, provided 
now that : 


|Z,,/Z2| <1, for a > 2 (79) 


Evidently the complete determination of the Z, can 
be carried out in this manner provided 


> |Z2] > ... > |Z, (80) 


It may be noted that the inequalities (80) imply the 
inequalities (46), thus guaranteeing the use of Eqs. 
(47). 


It is clear from Eqs. (70) that the degree to which the 
roots are separated in absolute value determines the 
rate of convergence of the process in a given prob- 
lem. 

In cases for which the convergence seems unsatisfac- 
tory it is convenient to revert to the determination of 
the stability equation itself. This can be accomplished 
through the use of m successive iterates of the type given 
by either Eqs. (59) or Eqs. (61). The procedure rests 
on the results that follow. 


Determination of Stability Equation by Iteration* 


Let it be imagined that the stability determinant 
given by Eq. (37) has been expanded; the resulting 
stability equation could be written in the form 


AZ 4+... + Ae =0 (81) 


* This is suggested in reference 11, p. 142, 


It is recalled that the Z, were defiried to be the roots of the stability equation; therefore, 


Za + + + A, = 0; 


aw1,2,...,% (82) 
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It may be stated that no assumption regarding the nature of the roots is intended in the development that 


follows. 


Let the first of Eqs. (82) be multiplied by a,a*'a,, the second by anza**a,™, etc. 


will be 


aya*'a,Z," + +... 


a,,a**a,Z," + a,a**a, + 


a,,a*"a, P a,,a*"a, Ay +. 
Adding the above equations, it is found that . 
a, + a, +. 


The resulting ‘equations 


+ A, + A, = 0 
A,-122 +a A, =0 


+ 


a, A,-1Z, + A, =0 


‘a, on 


Comparing the terms of this equation with Eqs. (67), it follows that 


Eqs. (84) evidently constitute a set of nonhomogeneous equations in the m unknowns Aj, Ao, ..., 


index 7 ranges from 1 to n. 


+ + A,a, =0 


(84) 


A,, since the 


In an analogous manner, the determination of the A; may be made to depend on the solution of the linear 


equations 
Am + + 


+ 4,-1.4, + 4,4, = 0 


(85) 


By matrix methods Eqs. (84) and (85) follow from the Cayley-Hamilton theorem, which states that 


+ + .. 


. + A,-.D + A, =0 


(86) 


As before, the A; are the coefficients of Eq. (81); D is the dynamical matrix whose elements are the D,,, while 


I is the identity matrix with elements 56,;. 
results obtained as Eqs. (67) and (83). 


It may be remarked that the theorem follows immediately from the 


Postmultiplication of Eq. (86) by an arbitrary column matrix gp with elements a; now yields Eqs. (84), while 
premultiplication of Eq. (86) by an arbitrary row matrix yp with elements A; results in Eqs. (85). 

An example illustrating the use of the modified iterative method just developed will be given in a later section. 
For the present another method of determining the stability equation is to be described. 


Classical Determination of the Stability Equation 


This process consists of the direct expansion of the 
stability determinant given by Eq. (37) to yield a 
polynomial in Z of degree n. The work involved can 
be conveniently organized through the use of Eq. (12). 
In the notation of that section, the D,; would correspond 
to the a1, say, and the (—4,,Z) would be taken for the 
by. The determinants corresponding to the right-hand 
member of Eq. (12) are then expanded according to 
Theorem 1 under the heading ‘‘Mathematical Prelimi- 
naries.”’ It may be noted that the coefficients of the vari- 
ous powers of Z are automatically separated by this 
procedure. The final result is of the form 


Diy — = (—1)"Z"* + + + 
(87) 


It may be noted for future reference that Eqs. (87) 
and (81) yield the result 
Ai = (88) 


The coefficients of the powers of Z which have been 
omitted from Eq. (87) are, of course, definite combina- 
tions of the Dy. Their general expressions are rather 


unwieldy and are therefore not given. In any par- 
ticular case, however, a complete literal expansion is 
highly desirable and may be carried through as indi- 
cated. 

For the purpose of numerical calculation, sets of de- 
tailed formulas leading progressively to the determina- 
tion of the A, for the stability equations of various de- 
grees should be prepared. It has been found that the 
utility of the classical procedure depends critically 
upon the extent to which the computing process is organ- 
ized. 

It will be noted that the last two sections have dealt 
only with the stability equation; the solution of Eqs. 
(36) requires in addition the consideration of two auxili- 
ary problems: 

(1) The solution of the algebraic equation of the 
nth degree with complex coefficients. 

(2) The solution of linear equations with complex 
coeficieuts. 

Each of these problems will be discussed in turn. 
For the present, however, a numerical example of the 
determination of the stability equation by the modified 
iterative process will be considered. 
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Example of Stability Equation Determined by Modified Iteration 
Consider the Equations 
(24 — Z)a,; + 10az + 11a; = 
7a, + (15 — Z)az — 17a3 = (1) 
23a, + 10a, + (12 — Z)as = 0 
Let A, be taken as (0, 0, 1); then 
10 11 
0} 7 0 9O 
1)|23 10 12 & 
23 10 12~ 
23)/24 10 11 552 230 253 
10}| 7 15 —17|| = || 70 150 —170 
12//23 10 12 276 120 144 , 
898 500 227 ~ A,® 
898}|24 10 11 21,552 8,980 9,878 f 
500) 7 15 —17|| = || 3,500 7,500 —8,500 
227/23 10 12 5,221 2,270 2,724 t 
tl 
30,273 18,750 4,102 ~ A, li 
Substituting these results into Eqs. (85), the following linear equations are obtained: 
30,273 + 898A; + 2342 + O0(A3) = 0 
18,750 + 500A; + + O0(A3) = 0 (2) 
4,102 + 227A, + 124.+ A; =0 
The advantage of using (0, 0, 1) for A; is now ap- 2, = 25; 2, = 25; Ze=l (5) 
— Sag It may be observed that, since Z; = Z:, the usual 
ti iterative determination of Z, would fail. 
It should perhaps be emphasized that the example Ea 
898A, + 234; = —30,273 (3) given illustrates the general case; i.e., in general, tip 
500A; + 104, = —18,750 A, would be taken to be (0, 0, ..., 0, 1) and A, firs 
However an additional simplification is possible. by (88). 
‘ : : procedure, only m — 2 linear equations need be solved 
It is recalled from Eq. (88) that in the general case A; = 4 , : 
to determine Ag, A3,..., A,-1; A, is then found from 
—D,; hence in the example under consideration it is the last of Eqs. (85) 
fount that It is highly recommended that the solution of the 
A, = —(24+ 15+ 12) = —51 n—2 linear equations be carried out by a matrix process 
Ths it follows from either of Bags. (3) above that to be described in the following section in all cases for 
—18,750 — 500A 
Ae = 675; ie., Ay = 0 = Solution of Linear Equations 
—18,750 + 25,500 The most important methods of solution are discussed 
10 = 675 in references 9, 10, 11, 12, 18, 15 and 16. The pro- 
cedures most useful for flutter purposes are believed 
The third of Eqs. (2) now yields A3; thus to be: 
As = —4,102 — 227 (—51) — 12(675) = —625 (1) Cramer’s rule as given by Eqs. (18). 
Hence, substituting into Bq. (81), the stability equation _ (2) | A matrix process described in reference 11, pp. [| *e 
for the system represented by Eqs. (1) is 125-130. 
Z? — 51Z? + 675Z — 625 = 0 (4) The solution by Cramer’s rule is convenient when lit- 
eral expressions for the unknowns in terms of the coef- 
the roots of this equation will be found to be ficients of the equations are desired. Such a solution 
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is highly desirable for the determination of amplitude 
ratios. The equations, it is recalled, are 
(Di 54jZ = 0; a= 2, 

The advantage of having detailed formulas expressing 
the a,, as functions of the D,, and the Z, is obvious. 
In this manner, the substitution of the values of Z, 
corresponding to a single set of D;; can be made the 
last step of the numerical computations. Moreover, 
many of the calculations involved are available if the 
stability equation has been obtained by the “‘classical’’ 
method. The formulas for calculating the a,;, should 
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3x, + 3x, = —7 
4x, + 2x2 + 5x3 — Vary = —5 


Let the coefficients of the unknowns be written in tabu- 
lar form; let the resulting array be denoted by M,; 


3 2 

3 1 -1 -3 

-1 —1 4 


Star any nonzero number in the first row of Mi; divide 
each number of the column containing the starred ele- 


therefore take advantage of this fact. 

The matrix process referred to is convenient and 
should be used whenever the numerical solution of a 
single set of equations is required. The process finds 


ment by the starred number and change the sign of 
each of the resulting quotients. Use the numbers ob- 
tained in this way as the first row of a new array, M;!; 
the remaining »—1 rows of M,' are to be the numbers 


its greatest application in conjunction with the modi- i's ug 
of Eqs. (84) or (85). It has been thought desirable, 1 2 1 
therefore, to include a numerital example illustrating 1 0 0 
the procedure in detail. Thus let the following set of Mi! = |/0 1 0 
0 0 1 


linear equations be considered: 


Calling the right members of the given equations /;, the following arrangement can be made to summarize the 
work up to this point: 


M, hy Mi} 
1 -I* 3 2 9 1 2 1 
3 i-l —@ -7 1 0 0 
4 2 5 —'Y/) —5 0 1 0 
~] 1. =} 4 3 0 0 1 


Each column of /;, except the column containing the starred element, is now used to construct a new column by mul- 
tiplying each of its elements into the corresponding row of M,! and vertically summing the results; thus using the 
first column of 1/;, it is found that 


1 0 oO -1 

Similarly it follows that 
2 1 21 2 2 4 2 

oO 1 0 oO -1 O 1 0 0 4 
-1 3%, 6 


The resulting new columns are now assembled to give a new array M;; it should be noted that the relative positions 
of the columns are preserved in going from M, to M2; thus it is found that 


4 2 -1 
11 3'/, 
0 2 6 


18 JOURNAL OF THE AERONAUTICAL SCIENCES—JANUARY, 1946 


The h, column is also used in the manner shown above to construct a new column, /2; thus 


2 1 
1 0 
3/0 O 1 


0 0 
0 
0 

2 13 12 


A nonzero element of the first row of M; is now starred; precisely the same process is now carried through with 


Mz as has been illustrated with M. 
Thus, beginning with the arrangement 


Ms In M; 
A 2 2 6 
6 11 13 
12 
new columns are constructed: 
6 14 24 2 6 7 12 3/2 6 7 12 
6|| 1 0|| = || 6 11 1 Of = 0}; 13 1 Oj = 0 
ie 0 2 0 1 0 2| 12 0 0 12 
20 «24 18 14 20 
As before it will be seen that the starred row of Mz now completes the second column 
of M; when this is done the fourth row of M will con- 
M; M;! tain a starred element; it is therefore to be completed 
20* 2 | | with zeros. Hence, it follows that 
24 1 4 
so that 0 0 
ll 3 2 
Continuing the pr in the er indicated, it is 
—7.6 found that 
and 1 4 20* 0 
M= —1* 0 0 0 


Consequently, M, is the single number —7.6; in this 
particular case /, is 0. Since the process that has been 
followed can go no further, it is clear that this phase of 
the work is complete. M/,, M2, M3, M, are now used to 
construct an “equivalent triangular matrix’”’ M, while 
hy, he, hz, hy are used to build a new column h. 

The starred row of M, now becomes the first column 
of M; the second row of M which now contains a 
starred element is to be completed with zeros; 
thus 


The new column h consists simply of the first elements 
of hy, he, hs, hy in that order; thus 


9 
2 
20 
0 


The original system of equations is known to be equiva- 
lent to the set 


x1 1 4 20* 0 9 
x O 0 0 _ 
Xs 3 2 24 —7.6*|| {20 
x 2 —1* 0 0 0 


Multiplying out the left member of this equation and 
setting the elements of the resulting column equal to 


zh with 


column 
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the corresponding elements of the right-hand column, 
it follows that 


— 3x3 + 2x, = 9 

4x; 2x3 = 2 

20x; + 24x3 = 20 

— 7.6x3 = 0 
Solving this ‘triangular’ set of equations it is found 
thatx3; = 0; x, = 1; x, = 2; x2 = —4. The answers 


_ may be checked in the original set of equations. 


Solution of the Stability Equation 


If the equation is a quadratic, the usual quadratic 
formula may be used. For the cubic and quartic equa- 
tions, the formal solutions of Tartaglia and Ferrari, 
respectively, have proved highly satisfactory. When 
dealing with complex coefficients, it is believed that the 
classical solutions referred to are definitely superior 
to the available methods of approximation. However, 
if the equation to be solved is of degree » greater than 4, 
it is necessary to determine »—4 roots by some method 
of approximation. The equation of degree » may 
then be reduced to a quartic after which Ferrari’s 
solution may be used. 

The computing procedures for the solutions referred 
to follow. 


SOLUTION OF THE GENERAL CuBIC EQuaTION 
AoZ* + + AZ + A; = 0 
Calculate in succession 
= '/3(Ai/ Ao); a2 = Ao); as = As/ Ao 
3/02; a1? — */oa2; ai(ai? — */2d2); 
P=, — = (a3/2) + ai(a;? — */202) 
U=-q+ VP*+ q’; change U to polar form 
w =U; change w back to rectangular form 
=—-—P/w; ww; wt; ww; wv 


where w = —(!/2) + j(V3/2) and = —(1/2) — 
i(V3/2). 


Roots are 
= ww + wy — ay 
Z3 = ww + w ay 
Checks: 


+ + = Ao 
ZiZ2Z3 = — do 
SOLUTION OF THE GENERAL QUARTIC EQuATION 
AoZ* + + + AZ + A = 0 
Calculate 


= "/4(Ai/ Ao); G2 =1/6(A2/ Ao); as = 1/4( As/ Ao) 


G37; + 3a.2; I = (a, + 3a2”) — 
+ a2* + as"; + a3” + 
J=+ (d204 + 2aia203) — + a3? + 
— a2; 2(a;? — a2); aide — 
T/12; (I/12)?; (1/12) 
J/8; (J/8)?; (5/8)? — (I/12)*; V(J/8)? — (I/12)* 
U = (J/8) + V(J/8)? — (1/12); change U to polar 
form 


w= U'*; change w to rectangular form 
(I/12)/w; + (I/12)/w 


Check: 
493 —-Ip —-J=0 
29; — 29; (a2 — 29)? 


(a2 — a; = — — 


= aig; ar; Zar 
v = 

Check: 
= (a — 29)*— ay Aat— am) +¢ 


+R = [2(a:? — a2) + ¢] = (v — 
V+R; V-R; 
Z, = (A — a) + V+R 


Checks: 


+ Z:Z3 + + Z2Z3 + + Z3Z4 = Ao 
+ + + = — Az/ Ao 
Z:Z2Z3Z, = do 


SUMMARY OF THE METHODS 


Three distinct calculation procedures for solving the 
flutter equations have been presented. A résumé of 
the manner in which each method has actually been 
used is now given. Table 1 is based on the use of the 
procedures in the manner to be described. 

It is understood that in all cases the given equations 
have already been reduced to canonical form. Thus 
an nth order stability determinant will always corre- 
spond to a stability equation of the mth degree. 

As a matter of general computing procedure it was 
thought desirable to have all calculations made by two 
computers working “‘in parallel.” The use of two people 
to do identical work is justified when the complete check 
that results and the greatly increased speed of each 
computer is taken into account. It may also be 
stated that the calculations were made by highly skilled 
personnel using late model Marchant and Friden cal- 
culators. In this connection it may be observed that 
from five to ten significant digits were carried in the 
computations. In making extended calculations, ac- 
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TABLE 1 
Time-Study Table 
Determination of the n 
Reduction Reduction Determination Sets of Amplitude Ratios 
of (n +1) of (m + 2) Determination of the of the n Corresponding to the n 
Number of Equations Equations nth Degree Stability Roots of the Roots of the Stability Calcula- 
Canonical ton ton Equatic — nth Degree Equation tion of 
Equations, Canonical Canonical Classical Modified Stability Classical Modified One 
n Equations Equations method iteration Equation method iteration Iterate 
3 0.50 1.25 0.70 0.70 2.50 1.00 1.00 0.25 
4 0.70 1.80 2.50 2.00 6.00 2.75 2.75 0.50 
5 1.00 2.50 12.50 4.50 8-10 7.50 7.50 0.80 
6 1.40 3.30 70 7.50 10-14 20 15 1.20 
7 1.80 4.25 450* 12.5 12-18 60* 28* 1.60 
8 2.40 5.30 3,500* 18.5 14-22 180* 48* 2.00 


* Estimated time; all numbers represent computing time in hours. 

Nore: It may be observed that the iterative process is represented only by the last column of data. The reason for this is that 
the number of iterates which must be computed to determine one root depends acutely on the separation of the roots and evidently 
cannot be stated in advance. In general, the process was found to be useful in conjunction with the classical method for n = 6— 
ie., one root was determined by iteration, thereby reducing the order of the determinant to be expanded from six to five. Forn < 
5 it was found that itt general there is no advantage to determining a root by iteration as a preliminary step in using the classical 


method. 


curacy of this order is necessary if well-defined results 
are to be obtained. It was felt that the uncertainties 
of original data are best reflected in the conclusions 
drawn from the final numerical results. 


A brief description of the use of each method now 


follows. 


Classical Method 


Literal expansions of the stability determinants of 
the third, fourth, fifth, and sixth orders were made. A 
thoroughly organized sequence of formulas building up 
to the coefficients of the associated stability equations 
was also prepared. Moreover, literal solutions for the 
amplitude ratios were made by Cramer’s rule. The 
complete set of formulas was then checked to avoid the 
duplication of any calculation. The computing work 
was organized for machine calculation. It is believed 
that a maximum of computing efficiency was achieved. 

The stability equations of the third, fourth, fifth, 
and sixth degrees were solved as previously indicated; 
for the quintic, one root was determined by Newton’s 
method of approximation, while for the sixth degree 
equation two roots were obtained by that process. 


Iterative Method 


This process was carried out as shown in the example 
given previously. The modification suggested there 
was used with the A, as well as with the a,”. 
However, the iterative process was never. used 
for the complete solution of a system of equations. 
For n = 4, one root was determined by iteration, after 
which the reduced set of three equations was handled 
by the classical method. For m = 5, one root was 
found by iteration, after which the classical procedure 
was followed. For m = 6 two roots were obtained by 
the iterative process. 


Modified Iterative Method 


This procedure was used in precisely the manner il-’ 


lustrated previously. Thus A, was taken to be 
(0, 0, ..., 0, 1), after which m successive iterates A,” 
were determined; A; was obtained from Eq. (88), 
after which the (m — 2) equations in As, As, ..., An—1 
were solved using the tabular process illustrated for 
solving linear equations, The last of Eqs. (85) was 
then used to determine A,. 

The resulting stability equations were solved as al- 
ready suggested. The amplitude ratios corresponding 
to the Z, were determined by the tabular process pre- 
viously referred to. P 


Comparison of Methods 


The basis for a just comparison of the methods out- 
lined would evidently be the total computing time re- 
quired for a ‘complete’ analysis. However, it will be 
readily admitted that the meaning of ‘‘complete” 
changes as an analysis proceeds. It has therefore 
been thought advisable to give the computing times as- 
sociated with each of the basic types of calculation. For 
a particular analysis the total times could be estimated 
from the data given. It should be noted that the figures 
given represent the computing times for each of two 
computers performing identical calculations independ- 
ently. 


Additional Considerations 


Before drawing any conclusions on the basis of the 
data presented certain other factors must be considered. 
In this connection it may be noted that the D,, of Eqs. 
(36) are functions of two auxiliary parameters—the air 
density p and the ‘“‘reduced wave length” 1/ko. Hence, 
at a definite altitude the flutter equations must be 
solved for a series of values of 1/ko sufficient to cover the 
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range of flying speeds of the airplane. However, even 
after this has been done, it is still desirable, on occasion, 
to solve the equations for values of 1/ko, intermediate 
to those for which solutions have been made. The ease 
with which this can be done is evidently an important 
consideration. 

A second point to be considered is that in general 

amplitude ratios are needed only in the following 
cases: 
(1) For zero air speed so that the choice of general- 
ized coordinates may be checked by comparing the cal- 
culated natural modes with ground vibration test 
data. 

(2) On rare occasions it is possible to check a calcu- 
lated flutter mode against flight vibration test data. 

(3) When a particular flutter mode seems to require 
a detailed study, the amplitude ratios are needed to pic- 
ture the mode. This sometimes makes it possible to 
carry out a more extended analysis of the mode with 
fewer degrees of freedom obtained by synthesizing those 
originally taken. 

A third consideration of paramount importance is 
the variation of the parameters defining the structure 
to be analyzed. The importance of flutter analysis in 
design is largely dependent on the extent to which the 
effect of such changes on the flutter characteristics of 
the structure can be readily determined. 

A fourth consideration is that of adding more degrees 
of freedom to an analysis already under way. It may 
be remarked that some indication of the reliability of 
an analysis can be obtained in this way. Thus the re- 
sults of a valid m degrees of freedom analysis would be 
substantially duplicated should n + 1 degrees of free- 
dom be used. 

Each of the three methods of computation will be 

reviewed relative to points one, two, three, and four in 
turn. 
(1) Solution of the Equations for Intermediate Values 
of 1/ko. If the stability equation has been found for 
several values of 1/ko, the real and imaginary parts of 
each of its coefficients may be plotted versus 1/ko, thus 
making it possible to interpolate the values of the 
coefficients for intermediate values of 1/ko. Evidently, 
the character of the graphs determines the accuracy of 
such interpolation. 

If an iterative solution has been made, the last iter- 
ate obtained for a given value of 1/k) would be used 
as the initial approximation for a neighboring value of 
1/ko. This process can be modified in obvious ways— 
for example, averaging iterates for values of 1/k) on 
either side of the desired value, etc. 

In actual practice, both the procedures described 
have been found to be satisfactory. There seems to 
be little reason to prefer one to the other. The com- 
puting times required have been found to be compa- 
table in most cases. 

(2) Calculation of Amplitude Ratios. This factor is 
not of too great importance, since, in general, the am- 


plitude ratios corresponding to a relatively small num- 
ber of roots are required. By virtue of this fact, how- 
ever, the superiority of the modified iterative method 
in this respect is greater than indicated by the time- 
study data, since the computing times given refer to the 
determination of m sets of amplitude ratios corre- 
sponding to the m roots of the mth degree stability equa- 
tion. 

(3) Variation of Parameters. In the classical pro- 
cedure it is clearly possible to isolate a parameter so 
that the coefficients of the stability equation become 
functions of that parameter. Moreover, several parame- 
ters can be partially isolated by a judicious expansion 
of the stability determinant. For example, in a recent 
flutter analysis the properties of the aileron were ef- 
fectively isolated so that changes in aileron ‘‘mass 
balance’’ and aileron “flapping frequency’ could be 
conveniently made. 

If either of the iterative processes is used, it is im- 
possible to isolate any parameter of the system. Each 
change must be investigated by first incorporating it 
into the original flutter equations. Next the canonical 
equations must be revised. If the iterative method is 
used, the last iterate for the original system of equations 
would now be used to suggest a first approximation for 
the modified system. 

If the modified iterative process is being used, a new 
set of linear equations for the coefficients of the stability 
polynomial would be calculated. The coefficients of 
the original stability equation would then be first ap- 
proximations for the solution of the revised set of linear 
equations. 

It is evident that if wide variations in any significant 
parameters are made, the iterative processes are not 
convenient. In this respect, the classical solution 
would be preferable, with the understanding, of course, 
that the parameters in question have been conveniently 
isolated. 

(4) Additional Degrees of Freedom. In the classical 
procedure the use of the last row or column in expand- 


ing the stability determinant for m + 1 degrees of free- 


dom makes it possible to utilize the calculations made 
in the m degrees of freedom expansion. It may also 
be noted that in such a case it is worth while to calcu- 
late first the amplitude ratios for the m degree case, 
since the work involved can also be used in the m + 1 
degree expansion. 

In the two iterative processes the results for the n 
degrees of freedom case would be used to suggest first 
approximations to the corresponding results for the 
n + 1 degrees of freedom problem. 

In practice this factor has not proved to be of great 
importance because of the fact that a sufficient number 
of properly chosen degrees of freedom seems to have 
been taken in all problems considered. Obviously, the 
fewer the degrees of freedom originally taken, the 
greater the possibility that additional degrees of freedom 
may bé later required. 


22 JOURNAL OF THE AERONAUTICAL SCIENCES —JANUARY, 1946 


CONCLUSIONS AND RECOMMENDATIONS 


On the basis of the comparison presented, the follow- 
ing conclusions are clearly in order: 

(1) The classical method and the modified iterative 
method are equally convenient when dealing with three 
canonical equations. It may be recalled that systems 
of four and five degrees of freedom often fall in this 

category when rigid-body motions are used as separate 
degrees of freedom. 

(2) The modified iterative method as described and 
illustrated in this paper is without question the best of 
the procedures given for solving sets of four, five, and 
six canonical equations. 

(3) If a set of four canonical equations is to be 
solved for several values of a parameter that can be 
conveniently isolated, it is advantageous to use the 
classical expansion process rather than the modified 
iterative method. This is the sole exception to con- 
clusion No. 2. 

It may also be noted that the determination of the 
stability equation by the modified iterative process does 
not involve prohibitive increases in computing time for 
n greater than six, However, the solution of the sta- 
bility equation does become considerably more time- 
consuming with increasing m. It is possible therefore 
that for some value of nm greater than eight the standard 
iterative method should be used as a preliminary step 
to determine a number of roots. The value of m in- 
volved cannot be stated until a thorough study of the 
approximate methods of solving algebraic equations is 
made. 

The author is aware of the fact that the conclusions 
stated are not those generally entertained by flutter 
analysts. The opinion has been commonly expressed 
that the iterative method should be used whenever the 
number of equations to be solved is greater than that 
for which the classical solution is feasible. The posi- 


tion held by the modified iterative method seems to 
have escaped attention. 
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A Simple Tabular Method of Calculating 


Deflections and Influence Coefficients of Beams 


N. O. MYKLESTAD* 
Consulting Engineer 


SUMMARY 


The present method of finding deflection curves and influence 
coefficients is easily applied to cantilever beams and to simple 
beams with or without overhangs. It is fast and accurate when 
a calculating machine is available, but slide rule accuracy is prob- 
ably sufficient for most practical purposes. The method possesses 
the additional advantage that the calculations can be performed 
in a routine manner by a computer who knows nothing about 
beam theory. 

The method assumes that the beam itself is weightless and car- 
ries a finite number of concentrated forces. 


INTRODUCTION 


6 be BASES OF THE PRESENT METHOD of analysis are 
the elastic coefficients of each section of the beam 
(between two consecutive forces), and these are defined 
as follows: Consider the section in question to be 
clamped at its left end so as to form a cantilever beam 
of length 1. When a unit force is applied to its right 
end, the slope and deflection at this end will be vp and 
dy, respectively, as shown in Fig. la. When a unit 
bending moment is applied to the right end of the 
clamped section, the slope and deflection at this end 
will be vy and dy, respectively, as shown in Fig. 1b. 
The quantities dy, dy, vr, and vy are called the “elastic 
coefficients.” 

The method consists of using the elastic coefficients 
to find the difference in the slopes at the two ends of 
the beam and also the difference in the deflections when 
the slope at one end is zero. The entire beam is then 
rotated about the support or end that has zero slope, 
until the boundary conditions at the other support are 
satisfied. 


DETERMINATION OF THE ELASTIC COEFFICIENTS 


For a uniform beam with flexural rigidity EJ, the 
values of the elastic coefficients are found from the 
simple beam theory to be 


= the.-? — (1) 
= dy = I?/2EI, Ibs.—! (2) 
dy = [/3EI, in. per lb. (3) 


If a section of the beam consists of two pieces, the 
elastic coefficients of which are both known, the elastic 
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coefficients of the two-piece section can be found from 
the formulas 


= + (4) 
Up = dy = vp’ + vp" + Ivy" (5) 
dy = dy’ + dp” + 1’ (2up” + (6) 


where the prime numbers refer to the right piece of 
the section and the double prime numbers refer to the 
left piece. 
If the bending stiffness curve varies smoothly alon 
the beam, it may be considered to be a straight line 
for each section, as shown in Fig. 2. For the right end 
of the section the cross-sectional moment of inertia is 
aly and for the left end, (a + 5)Jo, where 6 may be 
either positive or negative, but both a and a + 6 are 
always positive. Jo is a reference value for the moment 
of inertia and should preferably be taken approxi- 
mately equal to the smallest value of J for the entire 


4 


@) 


(6) 
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TABLE 1 
1 2 3 4 5 6 i 
(3) 2.3026(5) = b= 
n l a a+b (2) logio (4) In(4) (3) — (2) 
1 9.7 1.00 1.16 1.16 0.06446 0.1484 0.16 
2 12.3 1.16 1.56 1.345 0.12872 0.2964 0.40 
3 11.6 1.56 2.10 1.346 0.12905 0.2971 0.54 
4 15.1 2.10 2.70 1.286 0.10924 0.2515 0.60 
EI, = 8.7 X 108 
8 9 10 11 12 13 14 
vr X 106 = 
um X 10° = dy X 108 = dr X 10° = 
(1) (6)(8) (10)(11) ~ (2)(10) (8)(11)(13) 
n (7) 8.7 (7) — (2)(6) (8)(8) 8.7 8.7 
1 60.625 1.0341 0.01160 3,675.4 4.9005 pre 30.734 
2 30.750 1.0476 0.056176 945.56 6.1055 0.014836 49.582 
3 21.481 0.73356 0.076524 461.43 4.0587 0.026426 | 30.103 
4 25.167 0.72753 0.071850 633.38 5.2308 0.029115 53.345 
15 16 17 18 19 20 21 
S 1% 2[(12) (16) ant: X 10? = (14)(16) + (12)(17)] y X 10° = 
n F X 10- 108 for n — 1 + (9)(17)] (18), — (18) for n — 1 (20), — (20) 
1 0.150 0.150 | iia 0.73508 27.784 0 966.32 
2 0.250 0.400 1.455 4.7015 23.818 274.11 692.21 
3 0.300 0.700 6.375 12.219 16.300 595.79 370.53 
4 0.400 1.100 14.495 28.519 0 831.82 134.50 
b 31.105 28.519 966 . 32 0 
beam. The variable moment of inertia J for the sec- and the bending moment at station n is 
tion may then be written said 
M, = p> LS; (12) 
I = [a + (06/1) Wo (7) i=1 


and the elastic coefficients are 


In 
[s+ (24°) (10) 


a 


DETERMINATION OF DEFLECTIONS 


Cantilever Beam 

Consider first the case of a cantilever beam, as shown 
in Fig. 3. The positive directions chosen for forces, 
deflections, slopes, shear forces, and bending moments 
are such as to make all these values positive for the 
beam of Fig. 3. The shear force immediately to the 
left of the force F, is then 


(11) 


The nth section of the cantilever beam in its deflected 
position is shown in Fig. 4. From this figure the slope 
and deflection at station m + 1 can easily be obtained 
in terms of these same quantities at station m, giving 
the two equations 
(13) 
(14) 


= — UrnSn VunM,, 


When the slope and deflection at the free end of the 
cantilever beam are indicated by ¢ and 4, respectively, 
Egs. (13) and (14) give 
n—-1 
= — + vi Mi) 
+= 


n—1 


Ve = — + deiS; + dyiM;) 


(15) 
(16) 


where ¢ and 6 are determined by extending the summa- 
tions to the base of the cantilever beam and putting 
a = Y = 0 at the base. 


Example.—Consider the cantilever beam shown in 
Fig. 5, having a bending stiffness curve that may be 
considered to be trapezoidal for each section and with 
EI, = 8.7 X 10 Ibs.in.?, a, = 1.00, a2 = a + b = 
1.16, ds = @ + bp = 1.56, a, = a; + 0; = 2.10, and 
a4 + bs = 2.70. 
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TABLE 2 

1 2 3 4 5 6 7 8 9 10 11 

vr X 108 vp X 10° dr X 

= = 106 = 
du X 108 va X 10° da X 108 (1)'(6)" + 
(12) 10° ow X 108 = dr X 108 (1)'(5)" = = + 
l 2 1.5 ‘EI == (4)0) =(4)@) = (4)(8) (5)’ + (5)” 6)’ + (8) +7)" 

a’ 9.4 44.18 276.86 0.0426 0.40044 1.8821 11.794 


a” 10.1 51.005 343.43 0.0222 0.22422 1. 
419.90 0.0222 0.23976 1. 


10.8 58.32 


i 11.5 66.125 506.96 0.0128 0.14720 0. 
2’ 11.7 68.445 533.87 0.0128 0.14976 0. 


3.2400 0.62466 5.1221 60.518 
1323 7.6241 
2947 9.3218 

2.4362 0.38696 3.7309 51. 263 
84640 6.4891 
87610 6.8335 


1.1901 0.22368 2.0662 26.4712 


2” 8.8 38.72 227.16 0.00840 0.073920 0.32525 1.9081 

3 15 112.5 = 1,125.0 0.00840 0.12600 0.94500 9.4500 0.12600 0.94500 9.4500 
4’ 8.8 38.72 227.16 0.00840 0.073920 0.32525 1.9081 

4 ; 2.0498 0.22464 2.3751 30.943 
4” 9.6 46.08 294.91 0.0157 0.15072 0.72346 4.6301 

5’ 10.2 52.02 353.74 0.0157 0.16014 0.81671 5. 5537 

5 : 6.3040 0.57762 7.1207 104.09 

5” 9.8 48.02 313.73 0.0426 0.41748 2.0457 13.365 


The calculations for this case are all performed in 
Table 1, and the deflections themselves are obtained 
in Column 21 of this table. (20), in the operational 
index for Column 21 indicates the magnitude of the 
value in Column 20 for » = b—i.e., the last item in 


the column. 
Column 18 of Table 1 gives the value of the ex- 


n 
pression >> + and (18), must be equal to 
i=1 n 
since 0. The value of = ¢ — + 
vy.) is then obtained in Column 19. 


Simply Supported Beam 
Consider next the case of a simply supported beam, 


"as shown in Fig. 6. The positive directions chosen for 


forces, deflections, slopes, shear forces, and bending 
moments are now such as to make all these values posi- 
tive for the right side of the beam of Fig. 6. For the 
left side of this beam the slopes and shear forces will 
be negative, but the other quantities remain positive 
throughout the length of the beam. 
The shear force immediately to the left of the force 
F, is 
n 
(17) 
where R, is the reaction at the right end of the beam, 


considered positive upwards. The reaction R, at the 
left end of the beam is 


Ry = (1/L) (Fe (18) 


where L = >-J, is the total length of the beam and the 
summations without limits are to be extended over the 
entire beam. The reaction R, is then given by the 
expression 


n—1 


The bending moment at station m is 
M, = LS, (20) 
i=0 
The uth section of the simply supported beam in its 
deflected position is shown in Fig. 7. From this figure 
the slope and deflection-at station m + 1 are found as 


= A, — Urn Sn — UynM, (21) 
Ynt+1 = Vn + + drnSn + diunM, (22) 


The beam is now assumed to be rotated about the 
right support until the slope at this support is zero. For 
this rotated beam Eqs. (21) and (22) will give the slopes 
and deflections 


a,’ = — 2, + v4iM;) (23) 


n—1 ‘= 


where both a,’ and y,’ are always negative. In order to 

bring the left support back in place, the beam must now 

be rotated an amount ¢ = —y,’/L counterclockwise 

about the right support, where y,’ is obtained from Eq. . 
(24) by extending the summation over the entire beam. 

The slopes and deflections of the original beam are 


= a,’ (25) 
= ¢ + (26) 


Example.—Consider the shaft shown in Fig. 8, carry- 
ing five heavy discs of weights indicated in the figure. 
The shaft diameter varies in six steps, for which the 
flexural rigidities are 23.5 X 10° Ibs.-in.*, 45.0 X 10° 
Ibs.-in.?, 78.3 X Ibs.-in.*, 119.0 10° Ibs.-in.’, 
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TABLE 3 
2(3) = 2.325 2(2)(3) (10), 
1 2 3 4 5 7 s 9 10 1 
SX 10-* MX 10-3 2[(1)(8)—- X 10° 
Ro _ 2(1)(4) or X 10° = =[(4)(6) (5)(6)} (2) — 
n forn—1 FX 10-* 108 forn—1 dy X10* +(5)(7)] forn—1 (10) 
1.080 0 5.122 0.6247 5.532 60.52 0 0 
1 22.3 19.5 0.350 0.730 21.06 3.7309 0.3870 16.41 51.263 «42.53 392 
2 20.5 41.8 0.450 0.280 37.34 2.066 0.2237 25.34 26.47 202.5 «639 
3 15.0 62.3 0.550 0.270 43.08 0.9450 0.1260 30.51 9.450 727.4 661 
4 184 77.3 0.500 —0.770 30.03 2.375 0.2246 37.45 30.94 —1,147 576 
5 20.0 95.7 0.475 1.245 24.86 7.121 0.5776 42.94 104.1 366 
b 115.7 0 0 2,579 0 


63.7 X 10° Ibs.-in.?, and 23.5 X Ibs.-in.?, proceed- 
ing from right to left. The elastic coefficients for the 
sections of the shaft between discs, which sections are 
not the same as the constant diameter steps, are found 
in Columns 9, 10, and 11 of Table 2. Wherever a sec- 
tion of the shaft consists of two pieces of different diam- 
eters, the prime numbers refer to the right piece and 
the double prime numbers to the left piece. The elastic 
coefficients for the section may then be obtained from 
those of each piece by means of Eqs. (4) to (6). The 
deflections themselves are now found in Column 11 of 
Table 3. In this table Column 8 gives — a, +)’, accord- 
ing to Eq. (23), and Column 10 gives —y,’, according 
to Eqs. (23) and (24). 


Simply Supported Beam with Overhangs 


Assume now that the simply supported beam has a 
loaded overhang at each end. Each overhang may 
then be treated as a cantilever beam clamped at the 
nearest support. This will give the moments M, and 
M, at supports a and J, respectively, and the center 
span may then be treated exactly as a simple beam, 
except that the bending moment is decreased by an 
amount M, and the reaction R, is increased by the 
amount (M, — M,)/L. 

Eq. (17) remains the same as before, but Eqs. (19) 
and (20) now become 


n—-1 
R= (M, M,)/L + (1/L) (Fe (27) 
M, = M.+ (28) 


Eqs. (21) to (26) remain unchanged, and the tabular 
method of Table 3 can also be used for this case after 
performing the modifications indicated by Eqs.. (27) 
and (28). Thus the correct deflections are obtained for 
the center span. 

In order to obtain the correct deflections for the over- 
hangs, however, the right overhang has to be rotated 
an amount ¢ counterclockwise about the right support 
and the deflections thus obtained added to those found 
by treating this overhang as a cantilever beam. The 


left overhang must be rotated an amount —a’ — ¢ 
clockwise about the left support. The quantity — a’ 
is found at the bottom of Column 8 of modified Table 
3, and ¢ is given by (10),/(2), from the same table. 
The deflections for the entire beam are thus accurately 
determined. 


DETERMINATION OF INFLUENCE COEFFICIENTS 


Cantilever Beam 


Consider the cantilever beam of Fig. 3 but with all 
forces equal to zero. It is desired to find the deflection 
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TABLE 4 
1 2 3 4 5 6 7 8 9 10 
Z[(1)(4) + 
X anti’ X dr X 10® + GQnti”X Z[(1)(8)+ 

M’=2X(1) 10°+0% 108 = du X 108(2)] y’ X 108 = 106 = du X 10°] y” X 10% = 
n forn—1 1082)! (3), — (8) forn—1 °(5)»—(5) (7),—(7) forn—1 (9),— (9) 
1 9.7 0 4.9005 66.140 0 2,087.3 1.0341 2.5087 0 71.040 
2 12.3 9.7 21.168 49.873 672.29 1,415.0 2.0817 1.4611 29.235 41.805 
3 11.6 22.0 41.365 29.676 1,394.5 692.8 2.8153 0.7275 53.312 17.728 
4 15.1 33.6 71.041 0 1,858. 2 229.1 3.5428 0 65.810 5.230 
b 48.7 71.041 2,087.3 0 3.5428 71.040 0 


at station m due to a unit load at station m. This de- 
flection is denoted by 6,,, and from Maxwell’s Theorem 
it follows that dnn = 5am, Where dn, is the deflection at 
station m due to a unit load at station m. It is now 
further assumed that m < n. In Fig. 9a is shown the 
cantilever beam with a unit force at its free end, giving 
the slopes a’ and the deflections y’. In Fig. 9b is 
shown the same beam with a unit bending moment at 
its free end, giving slopes a” and deflections y”. In 
Fig. 9c the same beam is shown with a unit force at 
station m, and the problem now is to find the deflection 
5,m at station m, when m < n. Considering the part 
of the beam to the left of station m, it is easily seen 
that this part of the beam has the same loading as in 


_ Fig. 9d because the shear force and bending moment 


at station m are identical in the two cases and there is 
no load on this part of the beam. The deflection at 
station m due to the loading of Fig. 9d is easily seen 
to be 


m-—1 


5am = Yn’ — Yn” (29) 
Example.—Consider again the cantilever beam of 
Fig. 5, but with all the loads removed. It is now desired 
to find the influence coefficients for the four stations 
indicated. In Table 4 the values of y’ and y” are found 
by the method previously described. The values of 
the elastic coefficients are, of course, the same as 
were worked out in Table 1. In Table 5 the influence 
coefficients are found for the case when m < n. The 
values for m > n are found by Maxwell’s Theorem. 


Simply Supported Beam 


Consider now the beam of Fig. 6 but with all the 
loads equal to zero. It is desired to find the deflection 
at station m due to a unit force at station m; and, as 


TABLE 5 

m— 1 2 3 4 
0 9.7 12.3 11.6 
Lyn” 108 yn’ XK 108 = Yn’ — Yn”Zlm—1)108, S 2) 
1 71.040 
2 41.805 1,415.0 1,009.5 
3 17.728 692.8 520.84 302.78 
4 5.230 229.1 178.37 114.04 53.372 


Maxwell’s Theorem also holds for this case, it may be 
assumed that m < n. 

The simple beam is first made into a cantilever beam 
by removing its supports and clamping its left end. 
The values of y’ and y” are now found just as for the 
case of a cantilever beam. The right reaction of the 
simple beam due to a unit load at station m is 


m—1 


R, = 1 — (1/L) pe! 


where L is the total length of the beam and the right 
end section of the beam is now indicated by the sub- 
script 0 or a. It is also easily seen that the quantity 
m—1 


par’ is simply equal to the bending moment M,,’ of the 


cantilever beam with a unit force at its right end. The 
expression for R, may then be written 


Ra = 1 — (Mn '/L) (30) . 


Now load the cantilever beam with the force R, at 
its free end, as shown in Fig. 10a. The deflection at 
station m is then [1 — (M,,’/L)]y,’ upward, and the 
deflection of the right end is [1 — (M,,’/L) ]yq’ upward. 

Next, load the cantilever beam with a unit force at 
station m, as shown in Fig. 10b. The deflection at 
station m is now y,' — ¥y,”M,,’, and the deflection of 
the right end is y,,’.. When both these loads act simul- 
taneously on the cantilever beam, as shown in Fig. 10c, 
the upward deflection of station m will be B,M,,’, where 


2 = In” — /L (31) 
and the upward deflection of the right end is 
Am = [1 — (Mn'/L)]ye' — Ym’ (32) 


In order to make the right support level with the left, 
the entire beam is rotated clockwise about the left 
support an amount A,,/Z, which rotation will move 
station » down an amount [1 — (M,’/L)]Am. The 
total downward deflection of station m is then 


bam = [1 (M,’/L)]An — B,M,,’ (33) 


Example.—Consider again the shaft of Fig. 8, but 
without the discs. It is desired to find its influence co- 
efficients. In Table 6 the values of y’ and y” are ob- 
tained, using the values of the elastic coefficients ob- 
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TABLE 6 


1 2 3 4 5 6 
2[(1)(4) + 
M'= anti’X de X 108 
(1)... vy(2)} dy X10%2)) = 


n formn—1 X10* (3),—(8) forn—1 (5),—(5) 
a@ 19.5 0 5.122 113.62 0 10,107 
1 22.3 19.5 16.399 102.34 2,276.1 7,830.9 
2 20.5 41.8 27.816 90.924 4,682.3 5,424.7 
3 15.0 62.3 36.611 82.129 6,659.1 3,447.9 
4 18.4 77.3 56.347 62.393 7,959.3 2,147.7 
5 20.0 95.7 118.74 0 9,321.9 785.1 
b 115.7 118.74 10,107 0 
L = (2), = 115.7 
7 8 9 10 11 12 


(1)(8) AX 108 

xX10°= +dy y"X =(5)— BX 

— X 106 = (5)p 6 

(7) forn—1 (9-0) 

0.6247 1.5389 0 118.74 0 Riding 
1.0117 1.1519 35.131 83.609 572.68 15.926 
1.2354 0.9282 64.549 54.191 1,030.9 7.305 
0.8022 85.6438 33.097 1,216.9 3.297 
1.586 0.5776 98.621 20.119 1,206.8 1.556 


lor) 


2.1636 0 111.62 7.120 962.03 0.3343 
2.1636 118.74 0 0 0 
TABLE 7 
m—> 1 2 3 4 5 
Am X 10° 572.68 1,030.9 1,216.9 1,206.8 962.03 
M,,’ 19.5 41.8 62.3 77.3 95.7 
Xn 


1 (8mn — BaMn')108 (for m <n) 


ar 
8 
BS 
38 
8 
58 


0.33189 1.556 159.72 277.10 306.94 280.25 
53 182.77 134.30 


tained previously in Table 2. The values of A and B, 
given by Eqs. (32) and (31), respectively, are also 
worked out in Table 6. The value of y,’ is the quantity 
(5), of Table 6, and as y,’ = (5)y — (5),, from Column 
6, the value of A, may be written i 

= — (34) 
which is obtained in Column 11 of Table 6. The quan- 
tity B, is found in Column 12 of Table 6. In Table 7 
the influence coefficients themselves have been worked 
out for the case m < n. For brevity the quantity 
1 — (M,’/L) has been called x,/L, where x, is the dis- 
tance from the left support to station n. 


Simply Supported Beam with Overhangs 


For this case the subscripts R, C and L will be used - 


to indicate the right overhang, the center span, and 
the left overhang, respectively. 

Consider first the case where both stations m and n 
are in the center span. The influence coefficients for 


this case are determined as for a beam without over- 
hangs, which case was treated above. 

_ Consider next the case where station m is in the right 
overhang and station m either in the left overhang or in 
the center span. First, load the center span with a 
unit bending moment at its right support, as shown in 
Fig. lla. The value of Bg, in Fig. 1la is obviously 


Yen" — (Yen’/Le), given by Eq. (31). The value of | 


is 


= Bea/Le (35) 
the value of gz is 
= Aa" — (a¢a’ + Boa)/Le (36) 


and the deflection at station is Bggtcn/Le — Ben. 
It is now readily seen that when the unit force acts at 
station m of the right overhang, a distance xzm from 
the right support, the deflection at station mn of the 
center span is 

Sum = —Xprm(BoaXen/Le — Ben) (37) 


When station m is in the left overhang, a distance xz» 
from the left support, the deflection becomes 


Sum = — (38) 


Consider next the case where both m and n are in 
the right overhang, but m is closer to the right free end 
than is m. For this case the influence coefficient is 


Sam = + (39) 


where (dnm)z is the influence coefficient for the right 
overhang considered as a cantilever beam, clamped at 
support a. 

Next, consider the case where station m is in the left 
overhang and station m is in the center span. For this 
case the deflection of the center span due to a unit 
bending moment at the left support must first be ob- 
tained, as shown in Fig. 11b. . The deflection at station 
n for this case is obviously (Xen¥ea’/Le — Yen’)/Le, 
and the value of ¢,’ is 


= Yoa'/Le* (40) 


The deflection at station m of the center span due to a 
unit force at station m of the left overhang is then 


bam = — Yen’ \Xrm/Le (41) 


Finally, consider the case where both m and m are 
in the left overhang, but m is closer to the left free end 
than is m. The influence coefficient in this case is 


Sam = + (42) 


where (dnm)z is the influence coefficient for the left over- 
hang, considered as a cantilever beam clamped at 
support b. 

By means of Maxwell’s Theorem all influence co- 
efficients for the simply supported beam with over- 
hangs can now be found from the ones for each span 
separately, with the overhangs considered as cantilever 
beams, and Eqs. (37), (38), (39), (41), and (42). 
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General Instability of Semimonocoque 
Structures Under Bending 


TSUN KUEI WANG* 
Consolidated Vultee Aircraft Corporation 


SUMMARY 


Structural engineers are constantly faced with the difficulty 
of obtaining suitable formulas for the calculation of general 
instability of semimonocoque structures. The writer’ * * has 
written a series of papers in an effort to solve this problem. In 
this paper, a theory of buckling is presented in which an empirical 
correction is recommended. The theory takes into account the 
torsional rigidity of stringers and rings which other authors have 
not considered. The empirical correction introduced accounts 
for the effect of skin on the number of longitudinal buckling 
waves. Also, the formula as developed involves use of more 
of the geometrical characteristics and section properties of the 
structure than have the formulas developed by other authors. 
It is found that the theoretical predictions as modified by the 
correction factors and the experimental results are in agreement 
with regard to wave pattern and buckling load. 


PROPERTIES OF THE STRUCTURE 


E SEMIMONOCOQUE is the most commonly used 
| structure in present-day all-metal airplanes. This 
structure has three elements: stringers, transverse 
rings, and sheet covering. In this paper it is assumed 
that the stringers and rings are of constant cross 
section, the stringers are rigidly connected to the rings, 
and the sheet covering is securely fastened to both 
stringers and rings. Hereafter, the three elements are 
simply referred to as stringers, rings, and skin. 

In recent years the trend of development has been 
to increase the size of airplanes without increasing the 
dimensions of the three structural elements in the same 
proportion. With the type of structure that results 
from this development, general instability becomes 
vitally important. General instability is defined as 
the type of overall structural collapse in which stringers, 
rings, and effective skin buckle simultaneously, 


BrieF HistoRY OF PROBLEM AND DISCUSSION 
oF METHOD . 


Several analytical approaches have been made to this 
problem. Brazier, Dschou,® Taylor,’ Ryder,’ Hoff,® 
Heck,’ Timoshenko,” Nissen,'! and Cox!* have solved 
the problem with either one of the following two 
methods. One of them is the ‘‘equivalent shell” 
method (so designated at the California Institute of 
Technology"), and the other is the ‘‘equivalent truss’ 
method. In the former the flexural and torsional 
rigidities of the stringers and rings are assumed to be 
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distributed uniformly over the surface of the structure 
so that an orthotropic cylinder ftsults. In the latter 
the stringers, rings, and skin are considered as elements 
of a statically indeterminate truss. The normal load- 
resisting members of the truss are the stringers and 
rings, each with its proper portion of skin acting with 
it. The shear transferring members are the tension 
diagonals, each of which is formed by a suitable amount 
of skin in its own sheet panel. The values of the 
failing stress predicted by the various authors were 
checked by a comprehensive series of tests conducted 
at GALCIT."* It was found that the predictions 
were high by factors ranging as high as 12, except that 
those of Hoff were too low. The detailed discussion 
can be found in reference 14 and is not given here. 
There were also two late publications. One of them 
was Hoff’s theory,” which was obtained by a modi- 
fication of his original paper.* The other was the 
empirical formula of GALCIT. Both showed much 
better results. 

Obviously, testing each new design would be costly . 
and time-consuming. Therefore, a theoretical solution 
that involves suitable parameters and is in agreement 
with experiments is highly desirable. Since none of 
the present theories involves suitable parameters arid 
checks with experiment, further analytical research is 
necessary. As a whole the semimonocoque structure 
is anisotropic, and by nature the “equivalent truss” 
assumption is closer to the true picture. This would 
seem to render the “equivalent truss’ method pref- 
erable and is therefore used by the writer. 


ASSUMPTIONS 


In the present analysis the cross section of each 
stringer and ring is considered as augmented by an 
effective width of skin, and the buckling load of the 
structural grid made of effective stringers and rings is 
calculated. The additional restraint due to the tension 
diagonals'* of skin panels is taken into account later 
through an empirical correction. 

Further assumptions are as follows: 


(1) Stringers are identical in dimension and are 
equally spaced. The same is true of the rings. 

(2) End rings are absolutely rigid, while the 
stringers are hinged to them. 

(3) Rings fail in a plane normal to the axis of the 
cylinder. 
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(4) Stringers fail both radially and tangentially. 
(5) Hooke’s law is valid. ° 
(6) Neither local nor panel instability occurs. 


NOMENCLATURE 
The following nomenclature is used in the paper: 


(1) The geometrical characteristics group 
length of cylinder 
radius of cylinder 
number of intermediate rings 
number of stringers 
L/(p + 1) ring spacing 
2nrR/q stringer 
(2) The section properties group* 
E.J,, = flexural rigidity of stringers in the 
radial direction 
E,J,, = flexural rigidity of stringers in the 
tangential direction 


- cross-sectional area of stringers 
t thickness of skin 


(3) Buckling terminology group 
radial displacement of structure 


Ed; = flexural rigidity of rings (or frame) 
G,K, = torsional rigidity of stringers 
G,Ky = torsional rigidity of rings 


= 

w, = tangential displacement of structure 

Gm, = amplitude of buckling curve 

w, = effective width of skin 

m = number of half waves along stringer 
at buckling; this characterizes the 
wave length in the axial direction 

n = number of whole wave along ring at 


buckling; this characterizes the 
wave length in the circumferential 
direction 

(4) Loading nature group 


V = strain energy 

W = work done during buckling 

P, = applied force on the ith stringer 

Po = applied force acting on the most 
stressed stringer 

Poe = critical value of Po 

= w'E,I,,/L*, Euler load 

Pe = Po/P.y load ratio 


As shown in Fig. 1, the coordinate x is the distance 
in the longitudinal direction of the cylinder from its 
left end to right. The angle ¢ is measured from the 
vertical radius of the cross section. The sign con- 
vention used for these measurements may be seen from 
Fig. 1. 

It is assumed that the end rings of the structure do 


not change their shape during buckling. The de- — 


flections w, and w, along the end rings are therefore 
zero. 


* All the items include the effective width of skin except the 
last one. 


@, = 0 


when x = 0 and when x = L. 

The ends of the stringers are free to rotate at the end 
rings. Consequently, the end rings do not exert bend- 
ing moments on the stringers, and 


d*w,/dx? = d*w,/dx? = 


when x = 0 and when x = L. 

Experimental models used in the past had the 
stringers clamped to the end rings. However, in actual 
airplane structures perfectly rigid end restraint is 
hardly, if ever, realized. The assumption in the theory 
of simply supported ends is therefore conservative 
and, consequently, preferable. 

The general expression for radial deflections of the 
structure which satisfies the boundary conditions can 
be represented by a trigonometric series in which 

m= n=o 
w, = > Sin (mrx/L) cos nd (1) 
m=1,2,3,....n=2,3,4,.... 
where the number m starts from 2 because it represents 
an elliptical shape, the simplest type that can be 
obtained circumferentially by the structure during 
buckling. 

The tangential displacements are taken in such 
manner that the distortions of the rings are in- 
extensional, The discussion of inextensional deforma- 
tion of rings and shells was originated by Lord Ray- 
leigh'* and advanced by Timoshenko.” The condition 
for this requirement is 

w, = dw,/dd 


Thus, the tangential displacements of the structure 
under consideration can be expressed as: 


w, = sin n@ sin (mrx/L) 


m=1,2,;3,....m=2,3,4,.. 
(2) 


As long as the longitudinal stresses due to bending 
are below the yield point, the distribution of stresses 
is approximately according to the linear law. The 
external forces applied at the ends of the stringers may 
therefore be. assumed as 


P = Py cos 


For the calculation of the critical value of the force 
Po the Rayleigh-Ritz-Timoshenko energy method __ is 
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Fic. 1. General geometry of semimonocoque. 
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used. A buckling shape of the structure under con- 
sideration has been assumed as expressed in Eqs. (1) and 
(2). Corresponding to this shape, strain energy due to 
bending and torsion is stored in the rings, stringers, 
and skin. Moreover, at the same time work is done by 
the applied forces because the ends of the stringers 
undergo relative displacements during buckling. The 
critical value of Py is obtained by equating the total 
strain energy stored to the work done. The analytical 
expression is: 


(Visa + Visr) + (Visrz + Visr) = 
(Wise + Wis) (3) 


where the subscript 7 indicates the ith member, f refers 
to rings, s refers to stringers, r indicates in the radial 
direction, ¢ indicates in the tangential direction, B 
refers to strain energy caused by bending, and T refers 
to that caused by torsion. 

The terms in Eq. (3) can be calculated from the 


following integrals: 
dg 


(4) 
= So" 2, dd (5) 
Viers = (El er/2) So’ dx (6) 


View = (Esls:/2) So (O°w,/dx*)*, 4, dx (7) 


Vir = (G.K,/2R*) (d*w,/2¢dx)*, 4; dx (8) 
= (P,/2) (Ow,/ Ox) = 4; Ox (9) 
Wu = (P;/2) So (Ow,/ = 4; dx (10) 


EVALUATION OF INTEGRALS 


In the evaluation of the integrals use is made of the 
following notations: Let 


= diy sin aay sin ay sin 


L 

(11) 
Bin = Aig COS + cos Bang 

(12) 
Crt = Ama COS 24 + Ams COS 3G; + Gms COS 46; +... 

(13) 
Das = sin 36. + m sin 4¢,+ . 

(14) 


Emi = Sin 26; + Sin 3b, + 40m sin 46; + ... 
(15) 


Substitution of Eq. (11) in Eq. (1) yields: 


cos ng (16) 


The second partial derivative of Eq. (1) with respect 
to ¢ is: 


o¢? 


=- Sin cos nd 
L 


(17) 


The partial derivative of Eq. (1) with respect to x and 
¢ is: 


COS mn sin no (18) 
L 


With the aid of Eqs. (11) and (12), Eqs. (17) and 
(18) can be rewritten in the forms 


n=2,3,4,.. 


(19) 


= — («/L) nB,, sin nd (20) 


If Eqs. (16), (19), and (20) are substituted in Eqs. 
(4) and (5), it is found that the strain energies of bend- 
ing and torsion stored in the ith ring during buckling 
are 


In the same manner, by using Eqs. (6), (7), (8), 
(13), (14), and (15), the energies of bending and torsion 
stored in the ith stringer can be expressed as 


Views = (23) 
=1,2,3,... 

Var = (GK /ARL) 


Also, Eqs. (9), (10), (13), and (14) give the work 
done as: 


Wer = 


=1,2,3,.. 


(26) 


m? 


(27) 
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Substitution of Eqs. (21), (22), (23), (24), (25), (26), and (27) into Eq. (3) yields: 


E,l, i=p i=p n=o i=q m= 

Pass | 2— 1)24,,2 ‘C2 

Dini? | + 


m=o m =o 


2... Cal) +3 (cos¢; | (28) 
i=1 m=1,2,3,... i=1 m=1,2,3,... 


If the summation of the series is performed as shown in the appendix, Eq. (28) can be transformed into the 
following form: 


[#2 + ,L (n? + + 1)G,Ky > + 
m=1,2,3,...n=2,3,4,... RLq m=1,2,3,...0=2,3,4,... 
+ [1/n(n + }Gmn@mn+1 (29) 
m=1,2,3,...n=2,3,4,... 


By using the relation (p + 1)/g = (L/2mR)(c/b), where c/b is called panel aspect ratio, the Euler load P,, = 
rE, J,,/L?, the load ratio Pg = Po/P.,, and the nondimensional parameters 


mE] \b/\R4)’ WE wE,J,, |R? 


Eq. (29) can be simplified as below: 


m=1,2,3,...n=2,3,4, 
Pr = 


+ [1/n(m + 1) +1 


m=1,2,3,...n=2,3,4,... 


(30) 


CALCULATION OF CRITICAL VALUE OF Po 


To make Eq. (30) for Pg a minimum, a system of the and thus a system of linear equations with coefficients 


coefficients Gm, must be properly selected. This can @,, is obtained as follows: 


be done by taking the first derivative of Eq. (30) with 1 " 
respect to each coefficient and equating it to zero. Pp E + 41-2 + m? + 


Let 
2 
= (1) X + ont | + 
where [I] represents the numerator of Eq. (30) and l 
[II] the denominator. Then equation 0Pp/0dm, = 0 Pal + = 0 (33) 
gives: 


[L]/2a where = Ami = don = 0. 
2 = (31) It is seen that these homogeneous linear equations 
[IT]/Odmn iN Gm2, Gms, ... can be satisfied if all a,,, that correspond 


Eqs. (30) and (31) lead to the following result: to the unbuckled form of equilibrium of the cylinder 
are set equal to zero. To get nonzero solutions for 


[(: ») “4 dn, the constants Gms; - corresponding to a buckled 


shape of the cylinder, the determinant of the coef- 


1 r 1 . afi 1 ficients Gm, in Eq. (33) must be equal to zero. Thus 
(n n(n +1)]™"** an equation is obtained from which Pp,, can be calcu- 
(32) lated. 
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Eq. (33) may be written in the following manner: 


7 4+r_, 
0 0 g Pe + 
+ 4”) = 0 
m 
13 r 7 
0 0 -2(9+ = 0 
+ on) 
21 _5(16+r_, 13 
0 50 78 2( 16 “ft 75 Pe 0 = 0 
25H 167) 
m 
. 31 9 21 
39 Pe 2( 5 + 50 Pe 0 0 = 0 
+. 250) 


The more terms considered, the greater the accuracy of Pp. 


It was found that satisfactory approximation can 


be obtained from three equations corresponding to the assumption that only three coefficients—namely, @,,,, 


ANd Am4—are different from zero. 
the critical value of Pz is 


If the determinant of the coefficients of these equations is equated to zero, 


{[(4 + d)/4]m* + (9u/m*) + 4} {[(9 + d)/9]m* + (64u/m*) + 9v} 


(34) 


160{ [4 0)/4]m* + (9u/m*) + [(16 + d)/16]m? + (225u/m*) + 167) 
An approximation by taking four equations with four constants 


where m can be taken equal to 1, 2, 3, 4, .... 


Gms, ANA Ams Was also calculated but is not presented here. 


Comparison of the results of these calculations 


with the values given by Eq. (34) showed that the deviation was less than one-half of 1 per cent. 


It should be noted that Po is the force acting on two 
most stressed stringers, one being in compression and 
another in tension. For the compression stringer, the 
right-hand side of Eq. (34) shows in what proportion 
the stability of this stringer reinforced by rings is 
greater than the stability of an isolated bar. The 
magnitude of this proportional factor will vary when 
either one of the ratios I,,/Ig,, E_y/EDsr, 
K,G,/EJs, c/b, and L/R varies and also if m, the 
number of half waves into which the stringer buckles, 
varies. The number m (or the corresponding wave 
length L/2m) must be determined in such a manner 
as to make the right-hand member of Eq. (34) a mini- 
mum. This can be done graphically. The family of 
curves shown in Fig. 2 is a plot of Pz against » with 
m = 1, 2, 3, ... the parameter and with both ) and » 
assumed to be unity. For any given value of yu, the 
critical load corresponds to the value of Pz shown by 
the lowest of the curves. It is seen that for small 
values of u the curve m = 1 has the smallest ordinates. 
Hence, in this case the stringer buckles in one single 
half wave. Beginning with the point of intersection 
of curves m = 1 and m = 2, the second curve is the 
lowest and, consequently, the stringer buckles in two 


half waves. This holds true up to the point of inter- 
section of curves m = 2 and m = 3. The regions in 
which the stringer buckles in any other number of 
half waves can be defined in a similar manner. Figs. 
3 and 4 present only those portions of the curves that 
correspond to the minimum value of the ratio Pz with 
A’ = 1 and \ = 1.5, respectively, but the curves ar 
drawn for a number of values of v. ; 
The curves shown in Figs. 3 and 4 are correct only 
if the effect of the tension diagonals of the skin panels 


160 
140 
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Fic. 2. Pp, load ratio of buckling load P.. to Euler load Pex, 
against parameter uw with A = 1 and » = 1. ’ 
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on the number m of half waves of the stringers is not 
considered. Actually, the tension diagonals represent 
additional restraints that considerably affect the 
number m of half waves and, consequently, the buckling 
load, also. This effect is quite appreciable if the size 
of the skin panels is small. At present there is no way 
to evaluate this effect analytically, but the author was 
able to take it into account by means of correction 
factors established on the basis of the GALCIT tests. 
The correction factors are presented in Table 1 in 


TABLE 1 
Panel Size Relation Between m and m, 
bc = 2.5, in.? m, = 2m + 4.5 
be = 5, in.* m, = 2m + 3.5 
bc = 10, in.? m, = 2m + 2.5 
bc = 20, in.? me = 2m + 1.5 
bc = 40, in.? m, = 2m + 0.5 
bc = 80, in.? me = 2m — 0.5 


which m, indicates the corrected m. It is recom- 
mended that they be used in the following simple 
manner: 

Calculate first \, u, and vy and read from Fig. 3 or 
Fig. 4 the value of m that corresponds to the minimum 
buckling load. If \ = 1, Fig. 3 should be used. If 
\ = 1.5, Fig. 4 should be used. Determine m, with the 
aid of the formulas given in Table 1. Using the value 
of m, so obtained, together with the values of \, u, and 
v, find Pz with the aid of Figs. 5, 6, 7, or 8. (In these 
figures the subscript was omitted from m,.) In case 
\ is neither equal to 1 nor equal to 1.5, Figs. 3 and 4 
should be used. Simply follow the same procedure 
to obtain two P,’s, one of which corresponds to \ = 1 
and the other to X = 1.5. Then the value of Pp 
desired can be obtained by either interpolation or 
extrapolation. 
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Fic. 3. Pr, load ratio of buckling load P., to Euler load Pyy, 
against parameter » with \ = 1 and various values of ». 
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Fic. 4. Pr, load ratio of buckling load P., to Euler load P,., 
against parameter » with A = 1.5 and various values of ». 
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Fic. 5. Pr, load ratio of buckling load P., to Euler load Pru, 
against parameter » with X = 1 and various values of v and m. 
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Fic. 6. Pr, load ratio of buckling load P,, to Euler load Pex, 
against parameter u with \ = 1.5 and various values of vy and m. 
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Fic. 7. Pr, load ratio of buckling load P., to Euler load Px, 
against parameter » with \ = 1 and various values of vy and m. 


As a numerical example consider a cylinder for which 
1 u4=1,v = 1.4, and dc = 20 sq.in. From Fig. 
3, m is found to be 2. Then Table 1 gives m, = 5.5. 
Finally, with the aid of Fig. 5, the actual value of Pp 
is found to be 61. 


COMPARISON OF THEORETICAL WORK WITH 
EXPERIMENTAL DATA 


The comparison is based upon the comprehensive 
GALCIT tests on cylinders that failed by general 
instability. 

(1) Experimental results verify the prediction of the 
theory that the cylinder buckles according to a multi- 
lobe pattern in which the wave length is generally less 
than the total length of the cylinder. 

(2) In Table 2 experimental load ratios are com- 
pared to those predicted by the theory. It may be 
seen that there is good agreement between them. 


PROCEDURE OF APPLICATION OF THE THEORETICAL 
FORMULA TO ACTUAL PROBLEMS 


(1) Calculate the effective width 2w, of the skin 
acting with stringers by” 


2w,/c = Ore 


in which o,, = stress at the supported edge of the sheet 
panel and o,, = critical buckling stress of the plate. 

(2) Assume that the effective width of the sheet 
acting with the ring is equal to the total width of sheet 
between adjacent rings. 

(3) Calculate J;, J;,, and J,,.¥ 

(4) Calculate Ky, and K,.” 

(5) Calculate = 

(6) Calculate = 0.0103 
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Fic. 8. Pp, load ratio of buckling load P., to Euler load Pyy, 
against parameter uw with \ = 1.5 and various values of »v and m. 


TABLE 2 
Comparison with Experimental Data 


Model 

No. r m Me Prez, 
2 5.5 58.7 58 
 6§2.26..0.8 1.8 2 6.5 74 69 
27 #41.26 0.29 1.22 2 7.5 96 94.2 
0.92. 2 4.5 44 45.7 
St 0.37. 1.23 2 5.5 60 62.1 
32 1.25 0.58 1.31 2 6.5 79.5 83 
35 1.25 0.44 1.22 2 3.5 35.5 39.2 
0.7% 1.91 2 4.5 49 53 
7. 2 5.5 65 65.4 
88 1.25 0.91 1.49 2 7.5 101 93 
39 1.25 0.46 1.31 2 8.5 119 107 
41 1.24 0.11 0.98 2. 5.5 56.9 56.1 
42 1.24 0.19 1.00 2 6.5 73 66.5 
48 1.24 0.31 1.04 2 7.5 94 86.3 
45 1.24 0.22 1.00 2 4.5 42.2 45 
46 1.24 0.38 1.04 2 5.5 56 51.5 
47 1.24 0.62 1.12 2 6.5 77 77 
49 1.24 0.44 1.05 2 3.5 34 39 
5 1.2. 0.7. 1.2 2 4.5 46 46.2 
51 1.24 1.24 1.28 2 5.5 62.7 
54 1.25 2.26 1.13 3 8.5 123 122 
65 1.25 3.93 1.13 3 9.5 1651 155 
(7) Calculate » = 0.1013 [(KjG,/E,I,,)(c/b) + 

(K.G,/E,I.r) \(L/R)?. 


(8) If, in item 5, \ = 1, read value of m corre- 
sponding to wz and» from Fig. 3; if\ = 1.5, read m from 
Fig. 4. 

(9) Calculate panel size dc. 

(10) Obtain m, from Table 1. 
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(11) Enter the values of X, u, v, and m, (m, corre- (9’) Same as item 9. 
sponds to m in Figs. 5, 6, 7, and 8) into Figs. 5, 6, 7, or (10’) Obtain m, or m,’s from Table 1. 
8 to obtain Pz and then P,,. (11’) Enter the values of }’s, u, v, and m, (or m,’s) 


(12) The maximum ss aeeana buckling stress is into proper curves in Figs. 5, 6, 7, and 8 to obtain two 
Poer/ As. P's, one of which corresponds to \ = 1 and the other 


(8’) If, in item 5, \ is neither equal to 1 nor equal to \ = 1.5. Then an interpolation or extrapolation 
to 1.5, read value or values of m corresponding to » may be made according to the true value of ) in order 
and v from both Fig. 3 and Fig. 4. to obtain the desired value of Pp. 


Appendix 


For the purpose of simplifying the numerator and denominator in Eq. (28), the following steps are to be taken. 


It is evident that 
sin (mrx,/L) = sin [mxi/(p + 1)], cos (max,/L) = cos [mmi/(p + 1)] 


sin nd, = sin (2nmi/g), cos nd, = cos (2nri/q) 
In Eg. (28), the squares of Ain, Bin, Cry Daw and E,,, as shown in Eqs. (11), (12), (13), (14), and (15) contain 
finite summations of the following type: ay sin [mri/(p + 1)] sin [nri/(p + 1)]. 


An approximate mathematical solution i is y ar in order to render the final formula practical. It can be seen 
that for practical application, only a few terms of the general expression in Eqs. (1) and (2) need be considered. 
This makes both m and m small integers; while the structure under consideration is made up of a great number of 
stringers and rings so that integers p and g are relatively large. Therefore, it can be assumed that neither m nor 
n is greater than both (p + 1)/2 and g/2. Under this condition, the following conclusions can be obtained :?! 

(1) Either m + n or m — nis not a multiple of 2(p + 1), which causes 


[mxi/(p + 1)] sin + 1)] = 


when m # n. 
(2) Both m + nand m — mare not multiples of g, which causes 


j= 
sin (2mmi/q) sin (2nvi/q) = 
when m 
(3) Integer m is not a multiple of p + 1, which causes > sin’ [mni/(p + 1)] = (pb + 1)/2. 
i=1 


(4) Integer 2m is not a multiple of g, which causes }> sin? (2nri/q) = q/2. 


The above relations are also true with cosine products. Therefore, use of these relations leads to 


i=p n= m= n=o ) 


(n? — = + 1)/2] (n? — 1)dmn? 
i=1 n=2,3,4,... m=1,2,3,...=2,3,4,... 
t=p n=o m= 
n*Bin? = [e+ 1)/72] 
n=2,3,4,... 
i=q m=o m= n=@ 
= ‘mann? (35) 
t=1 m=1,2,3,... m=1,2,3,..."=2,3,4,... 
m'Dni? (q/2) (m4/n?)dmn* 
i=1 m=1,2,3,... m=1,2,3,...n=2,3,4,... 
m= m= n=o 
i=1 m=1,2,3,... m=1,2,3,...n=2,3,4,... 
Again the first term in the denominator in Eq. (28) can be simplified as follows: 
i=l m=1,2,3,... m=1,2,3,... i=l qd qd qd 
m?S (36) 
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. 


in which 
i=l q q q q qd q q 


Trigonometry gives the following identity: 


cos (2cri/g) cos (2mi/g) = */2 cos [2(¢ — 1)ri/g] + 1/2 cos [2(c + 1)ri/g] 


where c = 2, 3, 4, .... 


By use of this identity, Eq. (37) can be transformed in the following manner: 


i=q 


S= (ane cos + cos + aga cos + 


x cos + cos + ams cos + 
q q q 


q q q q 
+ + + + (ns + + + .. = (q/2) = +1 
Thus Eq. (36) becomes 
i=q@ m= m= 
DV (cosd mCn*) = (9/2) Onn +1 (38) 
i=1 m=1,2,3,... m=1,2,3,... n=2,3,4,... 
The second term in the denominator in Eq. (28) can be similarly simplified as follows: 
n=@ 
(cos > mi?) (q/2) [m?/n(n + 1) ]@mnmn +1 (39) 
i=1 m=1,2,3,... m=1,2,3,... #=2,3,4,... 


Substitution of Eqs. (35), (38), and (39) into Eq. (28) yields Eq. (29) as shown above. 
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Secant Modulus Method for Determining 
Plate Instability Above the 
Proportional Limit 


GEORGE GERARD* 
Republic Aviation Corporation 


SUMMARY 


The classical plate instability equation is considered in terms 
of strain rather than stress. Upon the assumption that the 
critical stress and critical strain are implicitly related by the 
stress-strain curve of the material, it is found that the secant 
modulus can be used in the plate buckling equation to predict 
critical stress above the proportional limit. 

The assumption is confirmed by tests on aluminum alloy Z 
and channel sections. 


SYMBOLS 
b = dimension of plate parallel to plane of loading, in. 
c = end fixity coefficient 
E = modulus of elasticity, Ibs. per sq.in. 
Ey = reduced modulus, lbs. per sq.in. 
Es = secant modulus, Ibs. per sq.in. 
E, = tangent modulus, lbs. per sq.in. 
k _ = plate instability coefficient 
K = — vy?) 
L = column length, in. 
DY = effective column length, L/+/c, in. 
t = plate thickness, in. 
x,y,2 = rectangular coordinates 
w = component of displacement in 2-direction, in. 
€ = strain, in. per in. 
o = stress, lbs. per sq.in. 
p = radius of gyration, in. 
v = Poisson’s ratio, 0.30 
T = ratio of E,/E 
= ratio of 
n = ratio of E;,/E 
Subscripts 


cr = critical 

pl = proportional limit 
f = flange 

w = web 


INTRODUCTION 


ere PRIMARY CONCERN in compression elements in 
aircraft structures usually centers about the in- 
stability characteristics of the element. These com- 
pression members are generally of such a configuration 
that they are subject either to primary column failure, 
or, if formed of thin sheet material, the sheet elements 


composing the cross section may be subject to local’ 


instability failure. 
Received December 29, 1944. Revised and resubmitted June 
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Below the proportional limit of a material, the prob- 
lem can be dealt with summarily by use of well-known 
stability equations. For primary instability, the Euler 
equation is used, 


o = crE/(L/p)? (1) 


and for local instability, the plate buckling equation 
as given by Timoshenko! is applicable: 


Ger = [kw?E/12(1 — v*)](t/b)? (2) 


Primary Instability Above Proportional Limit 


Above the proportional limit, the stability problem 
is inherently complex, since it is necessary to solve the 
equations for an anisotropic material. The theoretical 
solution for primary instability has been derived 
by von K4rman, who calculates a reduced modulus, 
E,, for a rectangular bar subject to bending stresses 
above the proportional limit. The reduced modulus 
as derived from the material stress-strain curve is 


E, = 4EE,/(E'* + E,*)? (3) 


Substitution of EZ, in place of E in the Euler equation 
(Eq. (1)) at the corresponding value of o will enable a 
plot of the theoretically correct column curve for the 
material. 

Sechler and Dunn? state that the column curves 
derived from the reduced modulus can be checked 
experimentally by specimens in which extreme care in 
manufacturing and testing has been exercised. How- 
ever, in practical structures testing, use of the tangent 
modulus, E,, in the Euler equation is in excellent 
agreement with test data for aluminum alloys. Fig. 1 


shows results of aluminum-alloy tests conducted at the - 


Aluminum Research Laboratories.*4 Thus, derived 
column curves based on the tangent modulus have 
real value in the practical primary instability problem. 


Plate Buckling Above Proportional Limit 


An applicable theoretical solution for buckling of thin 
plates above the proportional limit has not beén found, 
although attempted solutions have been published in 
the literature.» *°-? The local instability problem 
is of greater complexity than the primary instability 
case because of the stress condition in a buckled 
plate. 


38 


Fr 


Fic 


a 
A 
as 
stri 
Hor 
unl 
plat 
edg 
ing. 
buc 
deni 
T 
whe 
32 
in v 
benc 
in tl 
the | 
the 
modi 


DETERMINING PLATE INSTABILITY 39 


\ 

~“ \ \ © 75S-T REF4 
i — TANGENT MODULUS 

60 \ COLUMN CURVE 
3 
« 
3 
20 

° 10 20 30 40 60 60 70 60 90 100 


Fic. 1. Column curves showing agreement between tangent 
modulus and test data for several aluminum alloys. 


Y 
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Fic. 2. Deflection surface of: Jeft, column; right, plate supported 
at the unloaded edges. 


Assume a thin rectangular strip of sheet to be acting 
as a column. After the critical stress is reached, the 
strip will exhibit large deflections, these displacements 
occurring in the yz plane as shown in Fig. 2(left). 
However, the same strip, supported along one or both 
unloaded edges, will distort in both the xz and yz 
planes as shown in Fig. 2(right), since the unloaded 
edges are constrained to remain in the plane of load- 
ing. Thus, the fundamental difference between plate 
buckling and primary column failures is ay evi- 
dent. 

The differential equation of the deflected surface 
when stressed below the proportional limit is given by 
Timoshenko! as 


Et? & o'w 4 
+ 8 (a) 
in which the first term in the bracket accounts for 
bending in the yz plane and the third term for bending 
in the xz plane. The second term is concerned with 
the plate torsion which results from the interaction of 
the two bending components. 
Above the proportional limit, each of these terms is 
modified by a coefficient: 


dy? 
Et? 
where 7 is defined as 
= E,/E (5) | 


For the coefficient of the longitudinal bending term, 
Timoshenko! suggests the reduced modulus 


n=?Tt=E£,/E 
and, since the transverse bending is considered to be 
still in the elastic range 
1 
The evaluation of the coefficient for the torsion term 


is apparently difficult, but on the basis of a wide number 
of tests Kollbrunner® suggested using 


= (r+ V7)/2 


Using these values of the three coefficients, Lund- 
quist® arrived at an overall coefficient: 


= (r + 


Since r,is less than unity, it is evident that n' > + 
The only case in which an equality can exist is when 
E, = E, or below the proportional limit. 

Hence, using Lundquist’s notation for the effective 


modulus 7'£ above the proportional limit, the buckling . 


equation becomes 


While each term in Eq. (4a) has been modified for 
plasticity considerations, the solution of the k term has 
been accomplished by assuming a constant modulus 
of elasticity. Because of the complexities involved 
in obtaining a solution for k in which the modulus is 
variable, it is generally assumed that the k term does 
not change above the proportional limit. 

Because of certain assumptions made when deriving 
the coefficient for the torsion term, Lundquist suggests 
that 7»! be evaluated from the column curve of the 
material rather than from the reduced modulus. In 
the previous seetion, in which primary instability was 
discussed, it was shown (Fig. 1) that for aluminum 
alloys the tangent modulus column curves are in 
excellent agreement with test results. Consequently 
it is assumed that the column curve for the section 
shown in Fig. 2(/eft) is the tangent modulus column curve 
derived from the compression stress-strain curve. Thus 


E,/E 
and henceforth 
= (1! + 3 71)/4 


| 
| 
(6) 


Secant Modulus Method for Determining 
Plate Instability Above the 


Proportional Limit 
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SUMMARY 


The classical plate instability equation is considered in terms 
of strain rather than stress. Upon the assumption that the 
critical stress and critical strain are implicitly related by the 
stress-strain curve of the material, it is found that the secant 
modulus can be used in the plate buckling equation to predict 
critical stress above the proportional limit. 

The assumption is confirmed by tests on aluminum alloy Z 
and channel sections. 


SYMBOLS 


dimension of plate parallel to plane of loading, in. 
end fixity coefficient 

modulus of elasticity, Ibs. per sq.in. 

reduced modulus, Ibs. per sq.in. 

secant modulus, Ibs. per sq.in. 

tangent modulus, Ibs. per sq.in. 

plate instability coefficient 
kw?/12(1 — v?) 

column length, in. 
effective column length, L/ 
plate thickness, in. 
rectangular coordinates 
component of displacement in z-direction, in. 
strain, in. per in. 

stress, Ibs. per sq.in. 

radius of gyration, in. 

Poisson’s ratio, 0.30 

ratio of E,/E 

ratio of E,/E 

ratio of E;/E 

+ 3 -Vr1)/4 


in. 


432.3717 


Subscripts 
cr = critical 
pl = proportional limit 
f = flange 
w = web 


INTRODUCTION 


PRIMARY CONCERN in compression elements in 
aircraft structures usually centers about the in- 
stability characteristics of the element. These com- 
pression members are generally of such a configuration 
that they are subject either to primary column failure, 
or, if formed of thin sheet material, the sheet elements 
composing the cross section may be subject to local 
instability failure. 
Received December 29, 1944. Revised and resubmitted June 
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Below the proportional limit of a material, the prob- 
lem can be dealt with summarily by use of well-known 
stability equations. For primary instability, the Euler 
equation is used, 


o = crtE/(L/9)* (1) 


and for local instability, the plate buckling equation 
as given by Timoshenko! is applicable: 


= [kw?E/12(1 — v?)](¢/b)? (2) 


Primary Instability Above Proportional Limit 


Above the proportional limit, the stability problem 
is inherently complex, since it is necessary to solve the 
equations for an anisotropic material. The theoretical 
solution for primary instability has been derived 
by von Karman, who calculates a reduced modulus, 
E,, for a rectangular bar subject to bending stresses 
above the proportional limit. The reduced modulus 
as derived from the material stress-strain curve is 


E, = 4EE,/(E'? + E;’)? (3) 


Substitution of Z, in place of E in the Euler equation 
(Eq. (1)) at the corresponding value of o will enable a 
plot of the theoretically correct column curve for the 
material. 

Sechler and Dunn? state that the column curves 
derived from the reduced modulus can be checked 
experimentally by specimens in which extreme care in 
manufacturing and testing has been exercised. How- 
ever, in practical structures testing, use of the tangent 
modulus, £,, in the Euler equation is in excellent 
agreement with test data for aluminum alloys. Fig. 1 


shows results of aluminum-alloy tests conducted at the - 


Aluminum Research Laboratories.* 4 Thus, derived 
column curves based on the tangent modulus have 
real value in the practical primary instability problem. 


Plate Buckling Above Proportional Limit 


An applicable theoretical solution for buckling of thin 
plates above the proportional limit has not beén found, 


-although attempted ‘solutions have been published in 
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the literature.’ *°-’ The local instability problem 
is of greater complexity than the primary instability 
case because of the stress condition in a buckled 
plate. 
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Fic. 1. Column curves showing agreement between tangent 
modulus and test data for several aluminum alloys. 
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Fic. 2. Deflection surface of: Jeft, column; right, plate supported 
at the unloaded edges. 


Assume a thin rectangular strip of sheet to be acting 
as a column. After the critical stress is reached, the 
strip will exhibit large deflections, these displacements 
occurring in the yz plane as shown in Fig. 2(left). 
However, the same strip, supported along one or both 
unloaded edges, will distort in both the xz and yz 
planes as shown in Fig. 2(right), since the unloaded 
edges are constrained to remain in the plane of load- 
ing. Thus, the fundamental difference between plate 
buckling and primary column failures is clearly evi- 
dent. 

The differential equation of the deflected surface 
when stressed below the proportional limit is given by 


Timoshenko! as 


Et? 
- 12(1 — \Oxt - 
in which the first term in the bracket accounts for 
bending in the yz plane and the third term for bending 
in the xz plane. The second term is concerned with 
the plate torsion which results from the interaction of 
the two bending components. 

Above the proportional limit, each of these terms is 
modified by a coefficient: 


dy? 
Et? otw 
12(1 — (n + + % (4a) 


where 7 is defined as 
E,/E (5) 


For the coefficient of the longitudinal bending term, 
Timoshenko! suggests the reduced modulus 


n=rt=E£,/E 


and, since the transverse bending is considered to be 
still in the elastic range 


= 1 


The evaluation of the coefficient for the torsion term 
is apparently difficult, but on the basis of a wide number 
of tests Kollbrunner® suggested using 


= (r+ 


Using these values of the three coefficients, Lund- 
quist® arrived at an overall coefficient: 


nt = (r + 


Since zis less than unity, it is evident that »! > + 
The only case in which an equality can exist is when 
E, = E, or below the proportional limit. 

Hence, using Lundquist’s notation for the effective 
modulus 7!£ above the proportional limit, the buckling | 


equation becomes 
/(t\? 
™ 12(1 — (5) (6) 

While each term in Eq. (4a) has been modified for 
plasticity considerations, the solution of the k term has 
been accomplished by assuming a constant modulus 
of elasticity. Because of the complexities involved 
in obtaining a solution for k in which the modulus is 
variable, it is generally assumed that the k term does 
not change above the proportional limit. 

Because of certain assumptions made when deriving 
the coefficient for the torsion term, Lundquist suggests 
that 7! be evaluated from the column curve of the 
material rather than from the reduced modulus. In 
the previous seation, in which primary instability was 
discussed, it was shown (Fig. 1) that for aluminum 
alloys the tangent modulus column curves are in 
excellent agreement with test results. Consequently 
it is assumed that the column curve for the section 
shown in Fig. 2(/eft) is the tangent modulus column curve 
derived from the compression stress-strain curve. Thus 


r= 
and henceforth 
(11 + 38 
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Consequently, rewriting Eq. (6) so that the in- 
stability coefficient includes the numerical constants 
and Poisson ratio term, 


Ger/n' = KE(t/b)? (6a) 


The common procedure when using Eq. (6a) is to 
solve for ¢,,;/n! and then use curves such as shown in 
Fig. 3 to find the actual o,, as a function of ¢,,/n'. 
Sechler and Dunn? derive curves for various compres- 
sion yield values based on assumed parabolic short 
column curves. 
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Fic. 3. Theoretical variation of o¢- with o--/n! for 24S-T alumi- 
num alloy as given in reference 5. 


Equivalent Slenderness Ratio 


A rather simple solution of the buckling problem 
has been attempted by Langhaar*® and Hoff’ who em- 
ploy an “equivalent slenderness ratio’ concept. Using 
both stability equations, Eqs. (1) and (2), Langhaar 
assumes that the effective modulus for plates and 
columns is the same, and plots critical stress values 
versus L/mpv/c and b/t\/K as the abscissas. Hoff 
employs practically the same technique but uses the 
accepted column curves for the material. For rec- 
tangular sections, Fig. 2, this procedure is tantamount 
to assuming that 


Ger = E/12(1 — v*)](t/b)? (7) 


or 


= KE(t/b)? (7a) 


On the basis of the previous analysis of the differential 
equation, it is evident that the effecttve modulus for 
plate buckling is greater than that of columns. Conse- 
quently, the use of the equivalent slenderness ratio 
concept is questionable. 


CRITICAL STRAIN CONCEPT 


For convenience, the instability equation is written - 


as 
Oer/E = K(t/b)? (8) 
By applying the definition of strain in the elastic range 


e=a/E (9) 
Eq. (8) then becomes 
€cr = = K(t/b)* for Soy (10) 


Evidently the critical strain is solely a function of 
the geometry of the panel for a given material. While 
Poisson’s ratio enters into the evaluation of K, it has 
only a second order effect on this term. Consequently, 
panels of any structural material supported along the 
edges in the same manner and of identical geometry 
buckle at substantially the same strain. 

As previously noted, the modulus is no longer con- 
stant above the proportional limit and it is necessary to 
reconsider Eq. (2). In Fig. 4, the various moduli 
which are applicable above the proportional limit are 
defined on a typical stress-strain curve. The function 
which relates stress and strain directly is by definition 
the secant modulus; hence, Eq. (9) becomes in the 
yield region: 


e = 0/E, (9a) 


The assumption is now made that the critical stress 
and critical strain are implicitly related by the stress- 
strain curve of the material, i.e., the secant modulus 
relates the critical stress and critical strain. This 
assumption appears reasonable when it is considered 
that until instability initiates the element is essentially 
acting under block compression. 

Hence, based on this assumption, the plate buckling 
equation above the proportional limit can be written as 


= K(t/b)? 
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Fic. 4. Definition of the various moduli on a typical stress- 
strain curve for aluminum alloy. 
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and the restriction of Eq. (10) can be removed. In 
terms of critical stress, the expression of Eq. (9a) is 
used, 


Oe = KE,(t/b)? (11) 


Letting » = E,/E, the equation in coefficient form 
becomes 


Ser/n = KE(t/b)? (12) 


EXPERIMENTAL VERIFICATION 


To obtain experimental evidence of the validity of 
employing the secant modulus in the plate instability 
equation and substantiation of the asstumption that the 
critical stress and critical strain are implicitly related 
by the stress-strain curve of the material, the following 
test program was initiated. 


Specimens 

In the elastic range the theoretical buckling stress 
of plates simply supported along the unloaded edges 
can be determined from the basic plate buckling equa- 
tion. Experimental verification of the theoretical value 
is not often obtained because a simply supported edge 
condition is difficult to realize during a test. This is 
particularly true when the critical stress occurs in the 
plastic range. 

The use of Z or channel sections obviates the difficulty 
of supporting the test specimen, since the geometry of 
the section provides support for the flat plate elements 
comprising the cross section. In addition, the theo- 
retical buckling stress below the proportional limit for 


the Z or channel sections can be determined from the 


charts of reference 8, in which the theoretical work 
has been developed to a degree comparable with that 
of simply supported plates. Thus, use of these sections 
in the test program provides a convenient experimental 
method for duplicating boundary conditions specified 
by theory. 

Two sets of specimens were tested: one set was 
formed of 24S-O and the other of Alclad 75S-O. Both 
were subsequently heat-treated and aged to the -T 
temper. In this manner, increase in material proper- 
ties due to cold-working at the corners was avoided. 
The latter specimens were formed of Alclad material 
because bare 75S-T was not commercially available. 

Since a channel section and a Z section of the same 
dimensions have identical buckling stresses, each was 
included in the test program. The J, dimension (see 
Fig. 5) of all sections was nominally 2 in. To obtain 
different buckling stresses of the flange, the b, dimen- 
sion was varied in increments between 0.7 and 1.9 in. 
All specimens were designed to have the flange become 
unstable first. 


Test Method 


Each specimen was weighed and measured to ac- 
curately determine the cross-sectional area, thickness, 


abe 


STRAIN GAGES 


by 
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Fic. 5. Z and channel test specimens. 


and flange and web dimensions. These dimensions 
are listed in Table 1. On each specimen a slight 
variation in flange dimensions resulted from machining 
tolerances. The largest flange becomes unstable first, 
and thus two SR-4 strain gages were cemented back- 
to-back at the center of this leg for the instability 
investigation. 

The strain gages were used in determining both the 
critical stress and critical strain of the section in the 
manner shown in Fig. 6. The stress on the section as 
computed by load divided by area is plotted against 
the difference in the two strains existing on each side 
of the flange thickness. This difference in strain is a 
function of the bending stress of the flange which exists 
when buckling initiates. The buckling stress is de- 
termined as the point at which the bending strain 
becomes large with small increase in axial load. 

The critical strain is determined in a similar manner. 
The average of the two strain readings on each side of 
the flange thickness is the strain due to axial load. 
This value is plotted as the ordinate of Fig. 6. 

The specimens were 8 in. long to prevent premature 
failure by primary column action. Some specimens 
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TABLE 1 
Theoretical and Test Critical Strains and Stresses 
1 2 3 4 5 6 7 8 ) 10 11 12 
KwE(t/bw)?, 
Correction 


Kw, (t/bw)? KwE(t/bw)*, for Alclad, Test Test oc, 


Speci- by, buy t, 
men Material In In. In by/bw Ref. 8 X10-* Kips/In.2 Kips/In.2 Kips/In.? 
-1 24S-T 0.700 1.874 0.1285 0.374 3.58 4,706 180.3 ae Column failure 
-2 L 0.994 1.996 0.1275 0.498 2.64 4,083 115.3 i, 8,676 50.7 
-3 1.195 1.838 0.1280 0.650 1.75 4,844 90.7 ; 7,000 49.1 
-4 1.394 1.838 0.1275 0.758 1.38 4,816 70.6 5,080 46.0 
-5 1.597 1.963 0.1280 0.814 1.20 4,251 54.6 4,140 42.5 
-6 1.902 2.048 0.1275 0.929 0.95 3,869 39.3 arr 3,580 37.7 
-7 Ale. 75S-T 0.712 1.790 0.1280 0.398 3.40 5,112 166.9 144.1 Column failure 
-8 L 0.875 2.030 0.1265 0.431 3.11 3,881 115.9 100.2 7,606 72.5 
-9 0.987 2.050 0.1255 0.481 2.75 3,745 98.9 85.5 6,490 69.6 
-10 1.197 2.050 0.1255 0.584 2.08 3,745 74.8 64.7 7,000 66.5 
-11 1.403 1.998 0.1260 0.702 1.55 3,982 59.3 51.2 5,541 56.7 
-12 1.809 1.990 0.1260 0.909 1.00 4,007 38.5 33.3 3,370 40.0 
-13 24S-T 0.711 1.993 0.1280 0.357 3.68 4,122 162.3 7,304 52.5 
-14 os | 0.990 1.994 0.1275 0.496 2.65 4,083 115.8 9,018 50.0 
-15 1.205 2.014 0.1275 0.598 2.01 4,007 86.2 7,514 48.2 
-16 1.397 2.014 0.1275 0.694 1.57 4,007 67.3 5,750 44.0 
-17 1.591 2.030 0.1275 0.784 1.29 3,944 54.4 3,130 41.7 ’ 
-18 1.918 2.046 0.1280 0.937 0.93 3,881 38.6 say 3,400 34.7 ten 
-19 Ale. 75S-T 0.737 1.914 0.1275 0.385 3.50 4,436 149.1 128.9 14,950 74.7 are 
-20 | 0.866 1.932 0.1265 0.448 3.00 4,290 123.6 106.9 9,950 73.6 
-21 0.990 1.932 0.1265 0.512 2.53 4,290 104.2 99.1 8,400 71.5 
-22 1.213 1.908 0.1265 0.636 1.82 4,396 76.8 66.4 Not centered 
-23 1.436 1.916 0.1265 0.749 1.40 4,356 58.5 50.6 5,700 58.0 
-24 1.810 1.942 0.1265 0.932 0.95 4,238 38.7 | 33.4 4,030 40.8 
with short flanges appeared to have been influenced om wi 
by column action and were retested. Upon retest, the 
specimen length was reduced to 5 in. The Z and 5 fi ¥ belo: — 
channel sections were tested flat ended: a typical tf 
specimen at failing load is shown in Fig. 7. ty 3s Fi 
Critical Stress-Critical Strain Results as 
The test values of critical stress and critical strain om 2c 
are tabulated for each specimen in the last two columns , 004 32 
of Table 1. These test data are plotted on Fig. 8. - 
Included on this figure are the compression stress-strain . f 40 
curves of the Z and channel materials. A scatter 002 —— SPECIMEN- ~—— | 42 
band is drawn in the plastic range to depict the 24s-T C 43 
limits of variation in properties from the average R - 
curve. 40 47 
Within experimental accuracy, this figure confirms 50 
the basic assumption of the secant modulus method ” 
that the critical stress and critical strain are implicitly 30 
related by the stress-strain curve. For a given critical o-% rd 
stress, most of the test critical strain values are some- 2 ‘ 64 
what less than those predicted from the stress-strain gar 4 20 67 
curve. This is attributed to the experimental diffi- 69. 
culty involved in plastic strain measurements. The ng 
cement that bonds the strain gage to the stressed mem- to u “ 
ber has a tendency to creep at large strains. Hence, ° | bo 72. 
_ the critical strain as determined from strain gage €,- €, 74. 
measurements is effectively reduced from the true ", ’ a 


value. Fic. 6. Method of determining critical stress and critical strain. Sts 
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Fic. 7. Ling —— STRESS- STRAIN 
CURVE 


CORRELATION OF THEORIES WITH TEST DaTA 4 


The equations of each of the three methods for de- } I LL 
termining critical stress above the proportional limit 
are summarized: 

; Ser a 2 Fic. 8. Agreement between critical stress-critical strain test 
Lundquist, ni KE(t/b) (6a) data and stress-strain curve of the material. 
Langhaar-Hoff, a = KE(t/b)? (7a) Secant modulus, va = KE(t/b)? (12) 
T 
TABLE 2 
Calculation of Critical Stress by Various Theories 
1 2 3 4 5 6 7 . f 10 1l 
E,, @/®, E,, Fig. 8, 3/7, 
Fig.8, Fig. 8, Lbs. per Lbs.pr @/E @/E 3@% @+@ @ OO 
Kips/In.2 1078 Sq.In. X10® Sq.In. 108 Kips/In.2. Kips/In.2 Kips/In.? 
24S-T, E = 10.70 X 106 Ibs. per sq.in. 

29.7pl 2,780 10.70 10.70 1.000 1.000 3.000 1.000 29.7 29.7 29.7 

32.0 3,000 10.67 9.75 0.997 0.911 2.864 0.944 32.2 35.1 33.9 

36.0 3,500 10.29 6.36 0.961 0.594 2.312 0.727 37.5 60.6 49.5 

38.9 4,000 9.72 4.77 0.908 0.446 2.003 0.612 42.8 87.3 63.5 

40.9 4,500 9.09 3.44 0.850 0.322 1.702 0.506 48.1 127.2 80.9 

42.5 5,000 8.50 2.74 0.794 0.256 1.518 0.444 53.5 166.0 95.8 

43.7 5,500 7.94 2.16 0.742 0.202 1.348 0.388 58.9 216.4 112.8 

44.8 6,000 7.47 1.93 0.698- 0.180 1.273 0.363 64.2 248.3 123.3 

46.5 7,000 6.64 1.56 0.621 0.146 1.146 0.323 74.9 318.9 144.0 

47.9 8,000 5.99 1.30 0.560 0.122 1.048 0.292 85.6 394.2 163.8 

50.2 10,000 5.02 1.06 0.469 0.099 0.944 0.261 107.0 516.9 192.4 

53.0 13,000 4.08 0.80. 0.381 0.075 0.821 0.224 139.0 708.8 236.8 

: Alclad 75S-T, E = 9.60 X 108 Ibs. per sq.in. 

57.2pl 6,000 9.60 9.60 1.000 1.000 3.000 1.000 57.2 57.2 57.2 

61.2 6,500 9.42 7.37 0.982 0.768 2.629 0.849 62.3 79.7 72.1 

64.7 7,000 9.24 5.85 0.963 0.609 2.342 0.738 67.3 106.2 87.7 

67.2 7,500 8.96 4.21 0.933 0.439 1.987 0.606 72.1 153.3 110.8 

69.0 8,000 8.62 3.02 0.898 0.315 1.683 0.499 76.8 219.3 138.2 

70.4 8,500 8.28 2.26 0.863 0.235 1.455 0.423 81.6 299.1 166.6 

71.5 9,000 7.94 1.71 0.827 -0.178 1.266 0.361 ~ 86.5 401.5 198.1 

72.2 9,500 7.60 1.41 0.792 0.147 1.150 0.324 91.2 491.5 222.7 

72.8 10,000 7.28 1.24 . 0.758 0.129 1.079 0.302 96.0 563.5 241.1 

74.0 11,000 6.73 1.00 0.701 0.104 0.968 0.268 105.6 710.2 275.9 

75.0 12,000 6.25 0.91 0.652 0.095 0.926 0.255 115.0 787.8 293.9 

76.3 13,500 5.65 0.80 0.589 0.083 0.866 0.237 129.7 916.0 321.4 
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Each of these equations is identically equal to 
KE(t/b)?, although o,, as determined for each method 
is different since 7!, r!, and 7 are at variance above the 
proportional limit. The three coefficients have pre- 
viously been defined in terms of parameters derived 
from the compression stress-strain curve. From the 
average stress-strain curves of the Z and channel 
materials drawn on Fig. 8, and 
have been computed in Table 2 as a function of critical 
stress. 

In addition, the parameter K,E(t/b,)* has been 
calculated for each of the Z and channel sections tested 
in Table 1. In determining this parameter for Alclad 
material, it is necessary to account for the low yield 
strength cladding at the most highly stressed fibers. 
Tests by the Aluminum Company of America have 
indicated that buckling stresses are accurately pre- 
dicted if 93 per cent of the thickness is assumed effective 
in resisting instability. Consequently, the buckling 

eter was determined in this manner. 

For the specimens, the test o,, has been determined 
and thus the data of columns 9 and 12 of Table 1 and 
columns 1, 9, 10, and 11 of Table 2 permit an evaluation 
of the various methods for determining critical stresses. 
These data are plotted on Figs. 9 to show the agreement 
between test and theory. 

Excellent agreement between the secant modulus 
method and test data for the Z and channel section is 
obtained for 24S-T aluminum alloy. The agreement 
for Alclad 75S-T is not quite so good since the cladding 
introduces an indeterminacy that is difficult to account 
for. Although an effective thickness was used to 
compute K,,E(t/b,)?, the correction suggested by Alcoa 
is based largely on test data obtained in the elastic 
range whereas these tests are concerned principally 
with plastic strains. 


CONCLUSIONS 


The secant modulus method for determining critical 
stresses above the proportional limit has been con- 
firmed within experimental accuracy by tests on the 
aluminum alloys 24S-T and Alclad 75S-T. The agree- 
ment between theory and test data is good for strains 
up to 1.5 per cent, the maximum tested. It is reason- 
able to assuime that the secant modulus method is valid 
for all aluminum alloys in current use on aircraft struc- 

tures. 

To construct a chart of o,, as a function of KE(t/b)? 
above the proportional limit for a given material, it 
is merely necessary to multiply the abscissa. of the 
stress-strain diagram by the modulus of elasticity. By 
definition 


o/n = o/(E,/E) 
At any point on the stress-strain curve 
o/E, = 


SCIENCES—JANUARY, 1946 


2 
=: 
| 
20 
@ CHANNEL 
— THEORY 
10 
o 
° 20 40 60 60 100 120 140 


xipsyin? 


‘ Fic. 9a.. Agreement between test data and the various methods 
for 24S-T aluminum alloy. 
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Fic. 9b. Agreement between test data and the various methods 
for Alclad 75S-T aluminum alloy. 


and thus 
o/n = (13) 
An alternate procedure is to solve the equation 
€cr = K(t/b)? 
(Continued on page 48) 
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Stress-Strain Formulas 


WILLIAM R. OSGOOD* 
National Bureau of Standards 


ABSTRACT 


Most of the analytical expressions for stress-strain curves 
that have appeared in the literature and a few others are con- 
sidered critically. It is concluded that six have general value, 
and, of these, two are to be preferred in seeking a fit to stress- 
strain data. 


INTRODUCTION 


N VIEW OF THE INCREASING INTEREST in the analytical 
I representation of the stress-strain curve, it seems 
worth while to examine some of the expressions that 
have been proposed in the literature and to discuss 
briefly some of their advantages and disadvantages. 

First it may be well to consider some of the properties 
required and desired in any equation for a stress-strain 
curve. A prime requirement is simplicity, because an. 
equation is not to be set up merely for the beauty of its 
agreement with actual] stress-strain data. The equation 
is wanted for analysis of structural elements, beams, and 
columns, and any formula containing more than two 
parameters is almost sure to be either unnecessarily 
complicated or too complicated for practical use. On 
the other hand, one parameter is not enough to cover a 
wide range of materials with varying mechanical prop- 
erties and stress-strain curves of different shapes. Ele- 
mentary requirements are that the curve represented 
by the equation go through the origin and have a slope 
there equal to the modulus of elasticity. 

A formula should be as general as possible—that is, it 
should be adaptable to many different materials by 
varying the values of its parameters, thus precluding 
the necessity of trying several formulas to fit given data. 
It is desirable that the formula give physically possible 
values of the stress for all values of the strain and vice 
versa. It is also desirable that the formula be one that, 
when substituted in the integrals of plastic bénding 
and of column analysis, is not likely to lead to integrals 
that cannot be expressed in closed form. Finally, it is 
desirable that the parameters entering in the stress- 
strain formula be easily and accurately determinable 
for a good fit to experimental data. 

Single expressions covering both tension and com- 
pression will not be considered; absolute values of 
stress and of strain are assumed throughout this 
Paper. 

For strains up to about 1 per cent it is immaterial 
whether the stresses and the strains are based on the 
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original cross-sectional area and the original gage 
length or whether they are ‘‘true stresses” and ‘natural 
strains”,' but for larger strains the latter quantities 


should usually be used. 


POLYNOMIAL APPROXIMATIONS 


It is common in representing experimental data by 
empirical formulas to use polynomial approximations. 
Some of the empirical formulas that have been proposed 
to represent stress-strain curves are of this type, being 
special cases of the formulas 


Bo? + yo? + bot +... (1) 
o=e-+ be? + ce?+de4+... (2) 


where ¢ is the strain; o = s/E is the ratio of the stress, 
s, to the modulus of elasticity, Z; and 8, y, 6, ..., 0, 
c,d, ... are parameters to be determined for a good fit. 
By taking a sufficient number of terms, practically any 
stress-strain curve can be closely represented by either 
of these formulas. On the other hand, the stress-strain 
curves of many engineering materials can be adequately 
represented by one of these equations with the reten- 
tion of not more than two parameters. The question 
of which two to retain is not easily settled, however, 
and the very general expressions (1) and (2) may there- 
fore be dismissed as impractical. 


PROPOSED FORMULAS 


Many empirical formulas for representation of 
stress-strain curves have been proposed, and some of 
these will now be considered. Mehmke! has ksted a 
number of them. In the present discussion no formulas 
will be considered that represent curves that do not go 
through the origin or are discontinuous or have discon- 
tinuous slopes. 

The formulas containing only one parameter are 
in some instances special cases of more general 
formulas. Many of them were proposed for special 
materials, with no idea of generality in mind. It is 
therefore nut adverse criticism of the authors to point 
out that most of the formulas cannot be expected to 
apply generally. The quantities a, b, c, d, K, L, m, n, 
a, 8, and y, which occur in the formulas, are param- 
eters. 

t The true stress is the load divided by the actual cross-sec- 
tional area. The natural strain is the change in gage length 
divided by the actual gage length; it is the natural logarithm of 
one plus the strain based on the original gage length. 
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(a) The earliest relation expressed between stress 
and strain is Hooke’s law,? 


(3) 


This simple relation, a special case of Eq. (1) and famil- 
iar to everyone, hardly needs comment. It is accur- 
ate for low values of strain. Whether it is exact even 
for values approaching zero is at best doubtful. 

(b) In 1729 Bulffingeri® suggested 


e = Ko" (4) 


and this formula has been proposed several times since,’ 
in spite of the fact that the slope of the curve it repre- 
sents is not £ at the origin, unless K = nm = 1, in which 
case Hooke’s law, Eq. (3), results. Where large strains 
are involved, say of the order of 10 per cent, there may 
be some justification for its use. 


(c) Mehmke! attributes to Riccati‘ 


o = (5a) 
but as I read Ritcati, his considerations reduce to 
o = (1/a)(ew“t® — 1) (5b) 


in which e is the base of the natural system of log- 
arithms. Both formulas contain only one parameter, 
and the slope of the curve that the first one represents 
is zero at the origin. 


(d) Gerstner® proposed a formula of the form 
o=e-+ (6) 


a special case of Eq. (2). Since only one parameter 
enters, the formula cannot be expected to be generally 
applicable. Moreover, since ) must be negative, nega- 
tive values of o result when e > —1/b. 

(e) Poncelet® apparently was the first to suggest a 
formula that in its most general form, is largely free 
from objection. It may be written 


+ Be” — (7) 


It has the drawback that the parameters are difficult 
to determine for a best fit to stress-strain data. 


(f) Wertheim’ proposed 
e? = ao + Bo? (8) 
| The slope is zero at the origin, a serious limitation. 
(g) A formula proposed by Hodgkinson* may be 
written 
o =e-+ be? + ce? + det (9) 
It is a special case of Eq. (2), containing three param- 


eters, more than are desirable. Also, negative values 
of ¢ would result for certain ranges of values of e. 


(h) Formulas proposed by Cox® may be written 
e = a/(1 + ac) ’ (10) 


e=a+ + (11) 
o = e+ be? + (12) 


Eq. (10) contains only one parameter, and its use- 
fulness is therefore seriously limited. Eqs. (11) and (12) 
aré special cases of Eqs. (1) and (2). It is unlikely that 
either Eq. (11) or Eq. (12) will fit stress-strain data 
for many different materials, but in any given case it 
is not difficult to find out whether a good fit may be 
expected and, if so, to determine the values of the param- 
eters, as follows. For Eq. (11) plot (1/c)[(e/o) — 1] 
against ¢. If the result is approximately a straight line, 
Eq. (11) applies within the range of the observed values 
of ¢; and @ and y are the intercept and the slope of 
the line. 

For Eq. (12) plot (1/e)[(¢/e) — 1] against e. If 


the result is approximately a straight line, Eq. (12) 


applies within the range of the observed values of e; 
and 5 and ¢ are the intercept and the slope of the 
line. 

Eq. (11) will give values of e less than o, however, for 
a range of va ues of e if 8 or y is negative. Either } or c 
must be negative. If one of them is negative, Eq. (12) 
will give values of o greater than e for a range of values 
of e, and negative values of o for another range of values 
of e. If both } and ¢ are negative, Eq. (12) will give 
negative values of o for strains beyond a certain value 


of e. 
(i) Imbert” proposed a formula that suggests 
e = (1/a)(e*” — 1) (13) 
The formula contains only one parameter. 
(j) Formulas proposed by Hartig'' may be written 
og = (1/a)(e* — 1) (14) 
o = [e/(1 — e)jew (15) 
Each contains only one parameter—not enough to 
represent accurately stress-strain curves differing widely 
in shape. 
(k) Schiile!* proposed the formula 
o = Le” + be? (16) 
Here there are three parameters, but even so the slope 


is not E at the origin unless L = m = 1, in which case 
the formula reduces to Gerstner’s formula, Eq. (6). 


(1) Prager’s formula,'* 
o = ae + btanh [(1 — a)/ble (17) 


contains two parameters. It is of a form to be appli- 
cable to many materials where large strains are involved. 
It has the drawback that the parameters are not easily 
determinable for a best fit to experimental data (nor 
indeed to fit two points), and in problems of plastic 
bending even for simple cross sections it leads to inte- 
grals not expressible in closed form. 
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(m) Holmquist and Nadai'‘ proposed a double 
formula that, in its most general form, may be written 


e=a+ — a,)*, = apf (18) 
gp is intended to be the proportional limit. If, how- 


ever, g, is simply regarded as a third parameter, the 
formula becomes a powerful one. It could probably 
be made to represent accurately the stress-strain rela- 
tions of such materials as Alclad and some of the al- 
loys of magnesium that show a break in the lower part 
of the curve. Moreover, the formula would un- 
doubtedly give a more accurate description of a material 
that had been straightened by stretching than any other 
formula here presented. 

To obtain the parameters for the formula is some- 
what of a chore. A series of values of o, might be as- 
sumed and K and n then determined by plotting log 
(e — against log — Some values of would 
probably result in a straight-line plot for a reasonably 
large interval of e, and m and log K would be the slope 
and intercept of the straight line. If the corresponding 
value of a, led to a poor fit outside this interval, a com- 
promise between a,, n, and K would be necessary. 

The formula has the further disadvantage of being in 
two parts, thus adding to the complication of three 
parameters in using it for problems of plastic bending, 
for example. 


(n) Ramberg and Osgood” showed that the formula 
(19) 


is applicable to a wide variety of materials. This form- 
ula is a special case of Holmquist’s and Nadai’s formula, 
Eq. (18), from which it may be obtained by setting 
o, = 0. It is easy to find out whether Eq. (i9) 
is a good fit to stress-strain data and to determine the 
parameters if it is. If log (e — o) is plotted against log 
o and the result is a straight line, the formula will give 
a good fit, and the slope and the intercept of the line 
are n and log K. 

(o) A formula suggested by an unpublished for- 
mula by Rao and Leggett’® is 


e = 0+ £A(cosh ac — 1) 


This formula also has wide applicability, but the para- 
meters are not readily determined for a best fit. 

(p) Another formula’ suggested by a formula in 
Holmquist’s and Nadai’s paper" is 


e = ag + — 1) (21) 


This formula is of a form that probably makes it widely 
applicable where large strains are involved, but the 
parameters, again, are not easily determined for a best 
fit, even to fit two points. When a = 0, the formula 
reduces to Eq. (13). 

(q) A formula that may be found useful in cases 
involving large strains is 


= + be) 


(22) 
1+ ae 


- It is easy to find out whether the formula is a good fit 


to experimental data and to determine the parameters 
if it is. For a good fit a plot of (1/c) — (1/e) against 
e/o should turn out to be a straight line. If it does, a is 
the intercept on the axis of (1/0) — (1/e), and d is the - 
negative of the slope of the line. When 6 = 0, the 
formula reducés to Eq. (10). 


DISCUSSION 


The number of formulas that could be devised is of 
course unlimited. Some day empirical formulas may 
be replaced by a theoretical formula. Meanwhile, it 
should be possible to fit almost any stress-strain data 
adequately by means of one of the formulas in boldface 
type. Althvugh some formula with three parameters 
could be made to represent more stress-strain curves 
than any formula with two parameters, in general, the 
evaluation of three parameters for a best fit will be 
considerably more time-consuming than the evalua- 
tion of two, and a three-parameter formula will be 
much more cumbersome to apply in struciural analy- 
sis than a two-parameter formula. It does not seem 
worth while therefore to seek increased generality by 
means of a third parameter. Fig. 1 shows the curves 
represented by Eqs. (7), (19), and (20) with the param- 
eters in each case determined to make the curve pass 
through the points e = 0.006, « = 0.005 and e = 
0.0081, « = 0.0054. To the scale of the drawing the 
curves are indistinguishable from each other. Eqs. (17), 
(21), and (22) cannot be fitted to the two points chosen. 
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There is a significant difference between Eqs. (7), 
(19), and (20), on the one hand, and Eggs. (17), (21), 
and (22), on the other hand. As o increases for the 
former set, the slope of the o,e-curve approaches zero; 
whereas the o,e-curve represented by Eq. (17) ap- 
proaches the line 

o=aexb 


asymptotically, according as (1 — a)/b 2 0; the a,e- 
curve represented by Eq. (21) approaches the line 


e=ac— 8B 


asymptotically, unless (1 — a)/8 > 0, in which case 
the slope approaches zero; and the g,e-curve repre- 
sented by Eq. (22) approaches the line 


= (b/a)e + (1/a)[1 — (6/a)] 


asymptotically. These conditions may make Egs. 
(17), (21), and (22) more suitable than Eqs. (7), (19), 
and (20) when large (natural) strains and (true) stresses 
are important. 

It is suggested, when it is desired to represent a given 
stress-strain curve analytically, that Eqs. (19) and (22) 
be considered first, the former when small strains are 
important, the latter for large strains. Eq. (19) has 
been found to be widely applicable up to strains at 
least as great as 0.01. By plotting as suggested pre- 
viously, it is easy to find out whether these formulas 
may be expected to give good fits and to determine the 
parameters if they do give good fits. 
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‘Secant Modulus Method for Determining Plate Instability Above the 
Proportional Limit 


(Continued from page 44) 


and then use the stress-strain curve of the material to 
obtain the critical stress. 
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Fletcher Aviation Corporation 
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G & A Aircraft, Inc. 
General Aircraft Equipment, Inc. 
General Aviation Equipment Company, Inc. 
General Electric Company 
General Instrument Corporation 
General Motors Corporation 
AC Spark Plug Division 
Aeroproducts Division 
Buick Motor Division 
Cadillac Motor Car Division 
Chevrolet Motor Division 
Delco Products Division 
Delco-Remy Division 
Eastern Aircraft Division 
Fisher Body Division 
Frigidaire Division 
Harrison Radiator Division 
Research Laboratories Division 
Rochester Products Division 
The General Tire & Rubber Company 
The B. F. Goodrich Company 
The Goodyear Tire & Rubber Company 
The Gray Manufacturing Company 
Grumman Aircraft Engineering Corporation 
Guiberson Diesel Engine Company 
W. & L. E. Gurley 
Hanna Engineering Works 
Haskelite Manufacturing Corporation 
Hawaiian Airlines Limited 
Hayes Industries, Inc. 
The Hilliard Corporation 
The International Nickel Company 
Jack & Heintz, Inc. 
Jacobs Aircraft Engine Company 
Johns-Manville Sales Corporation 
Casey Jones School of Aeronautics, Inc. 
Kellett Aircraft Corporation 
Kenyon Instrument Company, Inc. 
Walter Kidde & Company, Inc. 


Kollsman Instrument Division, Square D Com- 


any 
Aircraft Corporation 
Langley Corporation 
Lavelle Aircraft Corporation 
Lear Incorporated 
Liberty Aircraft Products Corporation 
Lights, Incorporated 
Link Aviation Devices, Inc. 
The Liquidometer Corporation 
Lockheed Aircraft Corporation 
Longines-Wittnauer Watch Company, Inc. 
Magnaflux Corporation 
The Marquette Metal Products Company 
The Glenn L. Martin Company 
The W. L. Maxson Corporation 
Warren McArthur Corporation 
McDonnell Aircraft Corporation 
Menasco Manufacturing Company 
Micromatic Hone Corporation 
Minneapolis-Honeywell Regulator Company 
Moore Drop Forging Compan 
National City Bank of New York 
National Credit Office, Inc. 


The National Screw & Manufacturing Company 


The New York Air Brake Company 
Noorduyn Aviation Limited 
Bearings Corporation 
North American Aviation, Inc. 

North American Aviation, Inc., of Texas 


Northrop Aircraft, Inc. 

Northwest Airlines, Inc. 

Otto Aviation Corporation 
Owens-Corning Fiberglas Corporation 
Pan American Airways System 

The Parker Appliance Company 
Pennsylvania-Central Airlines Corporation 


Pesco Products Company Division, Borg-Warner 


Corporation 
Phillips Petroleum Company 
Pioneer Parachute Company, Inc. 
The Pure Oil Company 
Aviation Corporation 
J. P. Riddle Company 
John A. Roebling’s Sons Company 
Rohr Aircraft Corporation 
Romec Pump Company 
Roosevelt Field, Inc. 
The Ryan Aeronautical Company 
Sciaky Brothers 
Scott Aviation Corporation 
Shakespeare Products Company 
Shell Oil Company, Inc. 
Simmonds Aerocessories, Inc. 
Skydyne, Inc. 
Socony-Vacuum Oil Company 


Solar Aircraft Company 


Sperry Gyroscope Company, Inc. 
Square D Company 
Standard Oil Company of California 
Standard Oil Company (Indiana) 
Standard Oil Company of New Jersey 
Summerill Tubing Company 
Swedlow Aecroplastics 
Teleflex Limited 
The Texas Company 
Thompson Products, Inc. 
Tinnerman Products, Inc. 
Titeflex, Inc. 
Transcontinental & Western Air, Inc. 
Triplett & Barton, Inc. 
Union Carbide and Carbon Corporation 
Bakelite Corporation 
Haynes Stellite Company 
Linde Air Products Company 
National Carbon Company 
United Aircraft Corporation 
Chance Vought Aircraft Division 
Hamilton Standard Propellers Division 
Pratt & Whitney Aircraft Division 


Pratt & Whitney Aircraft Corporation of 


Missouri 
Sikorsky Aircraft Division 
United Air Lines, Inc. 
United States Aviation Underwriters, Inc. 
United States Rubber Company 
Dominion Rubber Company, Ltd. 
The Variety Aircraft Corporation 
Vickers, Inc. 
Victor Adding Machine Company 
Vidal Corporation 
The Waco Aircraft Company 
Warner Aircraft Corporation 
The Weatherhead Company 
Western Airlines, Inc. 
Westinghouse Electric Corporation 
Weston Electrical Instrument Corporation 
Wyman-Gordon Company 
Young Radiator Company 


SIXTH INTERNATIONAL CONGRESS 
FOR APPLIED MECHANICS 


Ata meeting initiated by Dr. J. C. Hunsaker and Dr. Th. von K4rman, acting as joint secretaries 
of the Fifth International Congress for Applied Mechanics, a committee was formed for the organiza- 
tion of the Sixth Congress. In line with the decision reached at Cambridge in 1938, it is proposed 
that the Sixth International Congress for Applied Mechanics be held in Paris, from September 22 to 
September 29, 1946. 

The invitations to the Congress are extended on behalf of 1’Académie des Sciences de 1’Institut 
de France, la Direction des Relations culturelles, le Centre national de la Recherche scientifique, 
l'Institut de Mécanique de la Faculté des Sciences de Paris, la Société francaise des Mécaniciens, 
and l’Association technique Maritime et Aéronautique. 


The Congress will meet at the Sorbonne. 
The Congress will be divided into the following Sections: 


I. Structures. Elasticity. Plasticity. 
II. Hydro- and Aerodynamics. Hydraulics. 
III. Solid Dynamics. Vibration and Sound. Friction and Lubrication. 
IV. Thermodynamics. Heat Transfer. Combustion. Fundamentais of Nuclear Energy. 
Besides the papers presented in these Sections, a number of General Lectures will be given on 
- subjects of current interest. The titles of these Lectures will be made known in a later notice. 


The French Organization Committee of the Congress is composed as follows: 


President: Henri VILLAT 
Secretary General: Maurice ROY 
Members: 


is de 

Paul DUMANOIS Joseph PERES 
Jules DRACH Ernest VESSIOT 


Those who desire to become members of the Congress are requested to inform the Secretary 
General as soon as possible of their intention to attend. They should also indicate at the same time 
whether they wish to present a paper, and in what Section. This is required in order to facilitate the 
preparation of the program and the issuing of further notices. 

Communications are to be addressed to the Secretary General of the Sixth International Con- 
gress for Applied Mechanics. 

Institut Henri-Poincaré, 11, rue Pierre-Curie, PARIS (V°) 


Checks or money orders can be made payable to the: 
Sixiéme Congrés International de Mécanique appliquée 
Institut Henri-Poincaré, 11, rue Pierre-Curie, PARIS (V°) 


It is recommended to make payments on the following account : Cheques postaux Paris N° 649 


(address given above). 

The authors who wish to present papers are requested to submit ‘he titles of their contributions 
as soon as possible. As at previous Congresses, contributions must be of scientific value and perti- 
nent to the program of the Congress. 

The authors are also asked to send to the Secretary General, before July 31, 1946, either the 
manuscript of their paper (maximum 5,000 words) or a summary of no more than 300 words. It is 
most desirable to have at least the summaries of all papers available to the members at the time of 
the Congress, in order to facilitate discussions. The Organization Committee hopes it will be able 
to have the summaries reproduced at the right time. However, it would be a great help if authors 
who have facilities to prepare a sufficient number of copies of their summaries, or manuscripts, would 
do so and send such copies to the Secretary or bring them to the meeting. 

The Proceedings of the Congress, containing the minutes of the sessions, the texts of the General 
Lectures, and the papers presented in the Sections, will be published with the least possible delay. 
The members of the Congress will have the right to _s these volumes by subscription. The surplus 
copies will be sold later on at a higher price. 

Adequate facilities will be provided for the members for their stay in Paris. Arrangements will 
be made known in due time. As at the previous Congresses, excursions and visits to various institu- 
tions, places of scientific, historic, and artistic interest will be arranged. 

The payment of an inscription fee of 300 francs is required from the members of the Congress; 
this fee can be paid immediately or not later than on the opening day of the Congress. Ladies ac- 
companying members will not be required to pay such fee, except those who wish to address the 


meetings or to submit papers. 


